sapply(c(
         "rjson", 
         "data.table", 
         "dplyr", 
         "ggplot2", 
         "stringr", 
         "purrr", 
         "foreach", 
         "doParallel", 
         "patchwork", 
         "lme4", 
         "lmerTest",
         "testit",
         "ggpubr",
         "latex2exp"
         ), 
       require, character=TRUE)
## Loading required package: rjson
## Loading required package: data.table
## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:data.table':
## 
##     between, first, last
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## Loading required package: ggplot2
## Loading required package: stringr
## Loading required package: purrr
## 
## Attaching package: 'purrr'
## The following object is masked from 'package:data.table':
## 
##     transpose
## Loading required package: foreach
## 
## Attaching package: 'foreach'
## The following objects are masked from 'package:purrr':
## 
##     accumulate, when
## Loading required package: doParallel
## Loading required package: iterators
## Loading required package: parallel
## Loading required package: patchwork
## Loading required package: lme4
## Loading required package: Matrix
## Loading required package: lmerTest
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
## Loading required package: testit
## Loading required package: ggpubr
## Loading required package: latex2exp
##      rjson data.table      dplyr    ggplot2    stringr      purrr    foreach 
##       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE 
## doParallel  patchwork       lme4   lmerTest     testit     ggpubr  latex2exp 
##       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE
sf <- function() sapply(paste0("./Functions/", list.files("./Functions/", recursive=TRUE)), source) # Source all fxs
sf()
##         ./Functions/General/FindDelay.R ./Functions/General/FindDelayAlt.R
## value   ?                               ?                                 
## visible FALSE                           FALSE                             
##         ./Functions/General/Utils.R ./Functions/Model/Models/RunMBayesLearner.R
## value   ?                           ?                                          
## visible FALSE                       FALSE                                      
##         ./Functions/Model/Models/RunMBayesLearnerDiffBeta.R
## value   ?                                                  
## visible FALSE                                              
##         ./Functions/Model/Models/RunMHybridDecayingQLearnerDiffBayes.R
## value   ?                                                             
## visible FALSE                                                         
##         ./Functions/Model/Models/RunMHybridDecayingQLearnerDiffBayesAlt.R
## value   ?                                                                
## visible FALSE                                                            
##         ./Functions/Model/Models/RunMHybridDiffDecayQLearnerBayes.R
## value   ?                                                          
## visible FALSE                                                      
##         ./Functions/Model/Models/RunMHybridDiffDecayQLearnerBayesAlt.R
## value   ?                                                             
## visible FALSE                                                         
##         ./Functions/Model/Models/RunMQLearner.R
## value   ?                                      
## visible FALSE                                  
##         ./Functions/Model/Models/RunMQLearnerDecayTo0Inits.R
## value   ?                                                   
## visible FALSE                                               
##         ./Functions/Model/Models/RunMQLearnerDecayToPessPrior.R
## value   ?                                                      
## visible FALSE                                                  
##         ./Functions/Model/Models/RunMQLearnerDecayToPessPriorDiffBeta.R
## value   ?                                                              
## visible FALSE                                                          
##         ./Functions/Model/Models/RunMQLearnerDecayToPessPriorDiffLR.R
## value   ?                                                            
## visible FALSE                                                        
##         ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPrior.R
## value   ?                                                          
## visible FALSE                                                      
##         ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPriorESAndEpsFixed.R
## value   ?                                                                       
## visible FALSE                                                                   
##         ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPriorNoCK.R
## value   ?                                                              
## visible FALSE                                                          
##         ./Functions/Model/ModelUtils/AltPlotSimTrainingCurves.R
## value   ?                                                      
## visible FALSE                                                  
##         ./Functions/Model/ModelUtils/aModelUtils.R
## value   ?                                         
## visible FALSE                                     
##         ./Functions/Model/ModelUtils/PlotAllWorstTestOptionsSimVsEmp.R
## value   ?                                                             
## visible FALSE                                                         
##         ./Functions/Model/ModelUtils/PlotSimEmpTest.R
## value   ?                                            
## visible FALSE                                        
##         ./Functions/Model/ModelUtils/PlotSimEmpTestNoSim.R
## value   ?                                                 
## visible FALSE                                             
##         ./Functions/Model/ModelUtils/PrepDfForModel.R
## value   ?                                            
## visible FALSE                                        
##         ./Functions/Model/ModelUtils/PrepDfForModelPROpt.R
## value   ?                                                 
## visible FALSE                                             
##         ./Functions/Model/ModelUtils/RecodeDfs.R
## value   ?                                       
## visible FALSE                                   
##         ./Functions/Model/ModelUtils/SimplePlotSimTrainingDelay.R
## value   ?                                                        
## visible FALSE                                                    
##         ./Functions/Model/OptFxs/RunOptTrainSIT.R
## value   ?                                        
## visible FALSE                                    
##         ./Functions/Model/OptFxs/RunOptTrainSITPR.R
## value   ?                                          
## visible FALSE                                      
##         ./Functions/Model/SimFxs/RunSimTrainSIT.R
## value   ?                                        
## visible FALSE                                    
##         ./Functions/Model/SimFxs/RunSimTrainSITForPR.R
## value   ?                                             
## visible FALSE                                         
##         ./Functions/Model/TrialWiseComps/aModelFunctions.R
## value   ?                                                 
## visible FALSE                                             
##         ./Functions/Model/TrialWiseComps/CalcBayesMean.R
## value   ?                                               
## visible FALSE                                           
##         ./Functions/Model/TrialWiseComps/CalcBayesStd.R
## value   ?                                              
## visible FALSE                                          
##         ./Functions/Model/TrialWiseComps/CalcBayesVar.R
## value   ?                                              
## visible FALSE                                          
##         ./Functions/Model/TrialWiseComps/DecayChoiceKernel.R
## value   ?                                                   
## visible FALSE                                               
##         ./Functions/Model/TrialWiseComps/DecayQVals.R
## value   ?                                            
## visible FALSE                                        
##         ./Functions/Model/TrialWiseComps/UpdateABMatrix.R
## value   ?                                                
## visible FALSE
DefPlotPars()
registerDoParallel(cores=round(detectCores()*2/3))

Load data from studies 1 and 2

s1_train <- read.csv("../data/cleaned_files/s1_train_with_delay_deident.csv")
  #read.csv("../data/cleaned_files/s1_train_deident.csv") %>% rename(ID=deident_ID) %>% relocate(ID)
s1_sit <- read.csv("../data/cleaned_files/s1_SIT_deident.csv") %>% rename(ID=deident_ID) %>% relocate(ID)
s2_train <- read.csv("../data/cleaned_files/s2_train_with_delay_deident.csv") #read.csv("../data/cleaned_files/s2_train_deident.csv") %>% rename(ID=deident_ID) %>% relocate(ID)
s2_sit <- read.csv("../data/cleaned_files/s2_sit_deident_corrected_names.csv") 

Harmonize variables and create some separate vars

# Study 2 harmonize  
s2_sit[s2_sit$condition=="cogn", "condition"] <- "cognitive" 
s2_train[s2_train$trial_within_condition <= 20, "block"] <- 1
s2_train[s2_train$trial_within_condition > 20, "block"] <- 2

s1_sit$probability <- factor(unlist(map(strsplit(as.character(s1_sit$valence_and_probability), "_"), 1)))
s1_sit$valence <- factor(unlist(map(strsplit(as.character(s1_sit$valence_and_probability), "_"), 2)))

s2_sit$probability <- factor(unlist(map(strsplit(as.character(s2_sit$valence_and_probability), "_"), 1)))
s2_sit$valence <- factor(unlist(map(strsplit(as.character(s2_sit$valence_and_probability), "_"), 2)))

Define paths and functions

# MLE results  
allp <- "../model_res/opts_mle_paper_final/all/"
# Empirical Bayes results 
bp <- "../model_res/opts_mle_paper_final/best/"
# Simulations (old folder name — includes emp bayes res)
sp <- "../model_res/sims_clean/sims_from_mle/"
# Read model function 
rm <- function(path, model_str) read.csv(paste0(path, model_str))

Behavioral Results for Learning and Test phase

Visualization

Summarize training mean and within-subject SEM

s1_train_summs <- s1_train %>% group_by(condition, ID) %>%  summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
s1_tr_summs <- Rmisc::summarySEwithin(s1_train_summs,
                        measurevar = "m",
                        withinvars = "condition", 
                        idvar = "ID")
## Automatically converting the following non-factors to factors: condition
s2_train_summs <- s2_train %>% group_by(condition, ID) %>%  summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
s2_tr_summs <- Rmisc::summarySEwithin(s2_train_summs,
                        measurevar = "m",
                        withinvars = "condition", 
                        idvar = "ID")
## Automatically converting the following non-factors to factors: condition
s2_tr_summs <- Rmisc::summarySEwithin(s2_train_summs,
                        measurevar = "m",
                        withinvars = "condition", 
                        idvar = "ID")
## Automatically converting the following non-factors to factors: condition
s1_test_summs <- s1_sit %>% group_by(condition, ID) %>%  summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
s1_te_summs <- Rmisc::summarySEwithin(s1_test_summs,
                        measurevar = "m",
                        withinvars = "condition", 
                        idvar = "ID")
## Automatically converting the following non-factors to factors: condition
s2_test_summs <- s2_sit %>% group_by(condition, ID) %>%  summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
s2_te_summs <- Rmisc::summarySEwithin(s2_test_summs,
                        measurevar = "m",
                        withinvars = "condition", 
                        idvar = "ID")
## Automatically converting the following non-factors to factors: condition
s2_te_summs <- Rmisc::summarySEwithin(s2_test_summs,
                        measurevar = "m",
                        withinvars = "condition", 
                        idvar = "ID")
## Automatically converting the following non-factors to factors: condition
s1_tr_summs$experiment <- 1
s1_tr_summs$phase <- "Learning"

s2_tr_summs$experiment <- 2
s2_tr_summs$phase <- "Learning"

s1_te_summs$experiment <- 1
s1_te_summs$phase <- "Test"

s2_te_summs$experiment <- 2
s2_te_summs$phase <- "Test"

comb_summs <- rbind(s1_tr_summs, s2_tr_summs)
comb_summs$experiment <- factor(comb_summs$experiment)

comb_te_summs <- rbind(s1_te_summs, s2_te_summs)
comb_te_summs$experiment <- factor(comb_te_summs$experiment)
a <- ggplot(comb_summs, 
      aes(x=experiment, y=m, group=condition, fill=condition)) + 
      geom_hline(yintercept=c(seq(.6, .8, .1)), size=.3) + 
      geom_hline(yintercept=.5, size=2) + 
      scale_fill_manual(values=c("purple", "orange")) +
      geom_errorbar(aes(x=experiment, ymin=m-se, ymax=m+se,
                        group=condition), 
                    inherit.aes=FALSE, size=1.5, width=.25, position = position_dodge(width = .3)) +
  geom_point(size=4, color="black", pch=21, position = position_dodge(width = .3)) + ga + ap + tol +
  xlab("") + ylab("Proportion correct") + ggtitle("Learning") + tp +
  ylim(c(.45, .82)) + theme(axis.text.x = element_text(size=30))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
a

b <- ggplot(comb_te_summs, 
      aes(x=experiment, y=m, group=condition, fill=condition)) + 
      geom_hline(yintercept=c(seq(.6, .8, .1)), size=.3) + 
      geom_hline(yintercept=.5, size=2) + 
      scale_fill_manual(values=c("purple", "orange")) +
      geom_errorbar(aes(x=experiment, ymin=m-se, ymax=m+se,
                        group=condition), 
                    inherit.aes=FALSE, size=1.5, width=.25, position = position_dodge(width = .3)) +
  geom_point(size=3, color="black", pch=21, position = position_dodge(width = .3)) + ga + ap + tol + tp + 
  #lp + Used to cut legend and then turned off so size is comparable 
  xlab("") + ylab("") + ggtitle("Test") +
  ylim(c(.45, .82)) + theme(axis.text.x = element_text(size=30))
b

perf_comb <- a + b
perf_comb

#ggsave("../paper/figs/fig_parts/perf_plot.png", perf_comb, height=6, width=10, dpi=700)

For resub, showing learning curves in Fig 2

s1_train_summs_lc <- s1_train %>% group_by(condition, ID, trial_within_condition, probability) %>%  summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition', 'ID',
## 'trial_within_condition'. You can override using the `.groups` argument.
s1_tr_summs_lc <- Rmisc::summarySEwithin(s1_train_summs_lc,
                        measurevar = "m",
                        withinvars = c("condition", "trial_within_condition", "probability"), 
                        idvar = "ID")
## Automatically converting the following non-factors to factors: condition, trial_within_condition, probability
s2_train_summs_lc <- s2_train %>% group_by(condition, ID, trial_within_condition, probability) %>%  summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition', 'ID',
## 'trial_within_condition'. You can override using the `.groups` argument.
s2_tr_summs_lc <- Rmisc::summarySEwithin(s2_train_summs_lc,
                        measurevar = "m",
                        withinvars = c("condition", "trial_within_condition", "probability"), 
                        idvar = "ID")
## Automatically converting the following non-factors to factors: condition, trial_within_condition, probability
s191_p <- ggplot(s1_tr_summs_lc %>% filter(probability=="90-10"), 
       aes(x=as.numeric(trial_within_condition), y=m, group=condition, color=condition)) + 
      geom_hline(yintercept = c(seq(.5, 1, .1))) +
      geom_hline(yintercept = .5, size=3) +
      geom_vline(xintercept = 20, size=1, color="gray57", linetype="dotted") + 
      geom_ribbon(aes(x=as.numeric(trial_within_condition), ymin=m-se, ymax=m+se), alpha=.3, color="gray57") + 
      geom_line(size=1, alpha=1) + 
      scale_color_manual(values=c("purple", "orange")) + ga + ap + 
  xlab("") + ylab("Proportion correct") + tol +
    ggtitle("90-10") + tp + theme(axis.text.x = element_blank(), axis.ticks.x=element_blank()) + 
  scale_y_continuous(breaks=c(.6, .8, 1)) + ylim(.48, 1.02)
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
s141_p <-ggplot(s1_tr_summs_lc %>% filter(probability=="40-10"), 
       aes(x=as.numeric(trial_within_condition), y=m, group=condition, color=condition)) + 
      geom_hline(yintercept = c(seq(.5, 1, .1))) +
      geom_hline(yintercept = .5, size=3) +
      geom_vline(xintercept = 20, size=1, color="gray57", linetype="dotted") + 
      geom_ribbon(aes(x=as.numeric(trial_within_condition), ymin=m-se, ymax=m+se), alpha=.3, color="gray57") + 
      geom_line(size=1, alpha=1) + 
      scale_color_manual(values=c("purple", "orange")) + ga + ap + 
  xlab("Stim iteration") + ylab("") + tol +
    ggtitle("40-10") + tp + 
  scale_y_continuous(breaks=c(.6, .8, 1)) + ylim(.48, 1.02)
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
s291_p <-ggplot(s2_tr_summs_lc %>% filter(probability=="90-10"), 
       aes(x=as.numeric(trial_within_condition), y=m, group=condition, color=condition)) + 
      geom_hline(yintercept = c(seq(.5, 1, .1))) +
      geom_hline(yintercept = .5, size=3) +
      geom_vline(xintercept = 20, size=1, color="gray57", linetype="dotted") + 
      geom_ribbon(aes(x=as.numeric(trial_within_condition), ymin=m-se, ymax=m+se), alpha=.3, color="gray57") + 
      geom_line(size=1, alpha=1) + 
      scale_color_manual(values=c("purple", "orange")) + ga + ap + 
  xlab("") + ylab("") + tol +
    ggtitle("90-10") + tp + theme(axis.text = element_blank(), axis.ticks=element_blank()) + 
   ylim(.48, 1.02)
s241_p <-ggplot(s2_tr_summs_lc %>% filter(probability=="40-10"), 
       aes(x=as.numeric(trial_within_condition), y=m, group=condition, color=condition)) + 
      geom_hline(yintercept = c(seq(.5, 1, .1))) +
      geom_hline(yintercept = .5, size=3) +
      geom_vline(xintercept = 20, size=1, color="gray57", linetype="dotted") + 
      geom_ribbon(aes(x=as.numeric(trial_within_condition), ymin=m-se, ymax=m+se), alpha=.3, color="gray57") + 
      geom_line(size=1, alpha=1) + 
      scale_color_manual(values=c("purple", "orange")) + ga + ap + 
  xlab("Stim iteration") + ylab("") + tol +
    ggtitle("40-10") + tp + theme(axis.text.y = element_blank(), axis.ticks.y=element_blank()) + 
   ylim(.48, 1.02)
s1_lcs <- s191_p / s141_p
s1_lcs
## Warning: Removed 1 row containing missing values (`geom_line()`).

s2_lcs <- s291_p / s241_p
s2_lcs

# ggsave("../paper/figs/fig_parts/s1_lcs.png", s1_lcs, height=6, width=5, dpi=700)
# ggsave("../paper/figs/fig_parts/s2_lcs.png", s2_lcs, height=6, width=5, dpi=700)

Statistical analyses

Create sum contrast codes

# Specify overt as baseline factor so that negative effect corresponds to worst performance  
s1_train$cond_cc <- factor(s1_train$condition, levels=c("overt", "cognitive"))
contrasts(s1_train$cond_cc) <- c(-.5, .5)
#head(s1_train$cond_cc)

s2_train$cond_cc <- factor(s2_train$condition, levels=c("overt", "cognitive"))
contrasts(s2_train$cond_cc) <- c(-.5, .5)

s1_sit$cond_cc <- factor(s1_sit$condition, levels=c("overt", "cognitive"))
contrasts(s1_sit$cond_cc) <- c(-.5, .5)

s2_sit$cond_cc <- factor(s2_sit$condition, levels=c("overt", "cognitive"))
contrasts(s2_sit$cond_cc) <- c(-.5, .5)

Regressions of condition differences in performance

Experiment 1 — Learning and test models

summary(m1_train_full_regr <- glmer(correct ~  
                                      cond_cc + scale(trial_within_condition) + 
                                      (cond_cc +  scale(trial_within_condition) |ID), 
          data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc + scale(trial_within_condition) + (cond_cc +  
##     scale(trial_within_condition) | ID)
##    Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  45054.2  45131.6 -22518.1  45036.2    39991 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -8.9157 -1.0686  0.4573  0.7047  1.2827 
## 
## Random effects:
##  Groups Name                          Variance Std.Dev. Corr       
##  ID     (Intercept)                   0.4466   0.6683              
##         cond_cc1                      0.4926   0.7018   -0.19      
##         scale(trial_within_condition) 0.1155   0.3399    0.82  0.01
## Number of obs: 40000, groups:  ID, 125
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.04191    0.06123  17.017  < 2e-16 ***
## cond_cc1                      -0.18054    0.06744  -2.677  0.00743 ** 
## scale(trial_within_condition)  0.36991    0.03300  11.209  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cnd_c1
## cond_cc1    -0.174       
## scl(trl_w_)  0.758  0.009
summary(m1_test_full_regr <- glmer(correct ~  cond_cc + (cond_cc|ID), 
          data=s1_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc + (cond_cc | ID)
##    Data: s1_sit
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   7854.4   7889.3  -3922.2   7844.4     7995 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.1729  0.1319  0.2789  0.5791  1.2316 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr 
##  ID     (Intercept) 1.558    1.248         
##         cond_cc1    1.972    1.404    -0.17
## Number of obs: 8000, groups:  ID, 125
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.5690     0.1193  13.150  < 2e-16 ***
## cond_cc1     -0.4569     0.1510  -3.027  0.00247 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## cond_cc1 -0.173

Experiment 2 — Learning and test models

summary(m2_train_full_regr <- glmer(correct ~  
                                      cond_cc + scale(trial_within_condition) + 
                                      (cond_cc +  scale(trial_within_condition) | ID), 
          data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc + scale(trial_within_condition) + (cond_cc +  
##     scale(trial_within_condition) | ID)
##    Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  48366.7  48445.0 -24174.4  48348.7    44151 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -10.0444  -1.0336   0.4405   0.6872   1.4792 
## 
## Random effects:
##  Groups Name                          Variance Std.Dev. Corr       
##  ID     (Intercept)                   0.6676   0.8171              
##         cond_cc1                      0.5648   0.7515   -0.07      
##         scale(trial_within_condition) 0.1676   0.4094    0.87 -0.04
## Number of obs: 44160, groups:  ID, 138
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.14912    0.07099  16.187   <2e-16 ***
## cond_cc1                      -0.17465    0.06851  -2.549   0.0108 *  
## scale(trial_within_condition)  0.39250    0.03732  10.517   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cnd_c1
## cond_cc1    -0.071       
## scl(trl_w_)  0.823 -0.034
summary(m2_test_full_regr <- glmer(correct ~  cond_cc + (cond_cc|ID), 
          data=s2_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc + (cond_cc | ID)
##    Data: s2_sit
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   8426.3   8461.7  -4208.1   8416.3     8827 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.2283  0.1250  0.2301  0.5751  1.5541 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr 
##  ID     (Intercept) 1.850    1.360         
##         cond_cc1    2.467    1.571    -0.15
## Number of obs: 8832, groups:  ID, 138
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.6720     0.1235  13.543   <2e-16 ***
## cond_cc1     -0.3543     0.1586  -2.234   0.0255 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## cond_cc1 -0.156
car::vif(m1_train_full_regr)
##                       cond_cc scale(trial_within_condition) 
##                      1.000073                      1.000073
car::vif(m2_train_full_regr)
##                       cond_cc scale(trial_within_condition) 
##                      1.001173                      1.001173

Robustness check 1: Same models with just random intercepts instead of both random intercepts and slopes

Experiment 1

summary(m1_train_ri_only_regr <- glmer(correct ~  
                                      cond_cc + scale(trial_within_condition) + 
                                      (1 |ID), 
          data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc + scale(trial_within_condition) + (1 | ID)
##    Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  46071.1  46105.5 -23031.6  46063.1    39996 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.6755 -1.1011  0.5027  0.6818  1.2499 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  ID     (Intercept) 0.322    0.5675  
## Number of obs: 40000, groups:  ID, 125
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    0.95856    0.05213  18.388  < 2e-16 ***
## cond_cc1                      -0.14312    0.02271  -6.302 2.93e-10 ***
## scale(trial_within_condition)  0.28711    0.01147  25.021  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cnd_c1
## cond_cc1    -0.006       
## scl(trl_w_)  0.026 -0.008
summary(m1_test_ri_only_regr <- glmer(correct ~  cond_cc + (1 |ID), 
          data=s1_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc + (1 | ID)
##    Data: s1_sit
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   8080.6   8101.5  -4037.3   8074.6     7997 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4964  0.1516  0.3789  0.6174  1.2182 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  ID     (Intercept) 1.299    1.14    
## Number of obs: 8000, groups:  ID, 125
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.40802    0.10812  13.023  < 2e-16 ***
## cond_cc1    -0.31767    0.05567  -5.707 1.15e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## cond_cc1 -0.022

Experiment 2

summary(m2_train_ri_only_regr <- glmer(correct ~  
                                      cond_cc + scale(trial_within_condition) + 
                                      (1  | ID), 
          data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc + scale(trial_within_condition) + (1 | ID)
##    Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  49635.6  49670.3 -24813.8  49627.6    44156 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.4739 -1.0679  0.4853  0.6668  1.2916 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  ID     (Intercept) 0.4132   0.6428  
## Number of obs: 44160, groups:  ID, 138
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.03285    0.05597  18.452  < 2e-16 ***
## cond_cc1                      -0.15224    0.02195  -6.937 4.02e-12 ***
## scale(trial_within_condition)  0.27618    0.01108  24.921  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cnd_c1
## cond_cc1    -0.007       
## scl(trl_w_)  0.024 -0.008
summary(m2_test_ri_only_regr <- glmer(correct ~  cond_cc + (1 |ID), 
          data=s2_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc + (1 | ID)
##    Data: s2_sit
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   8769.8   8791.0  -4381.9   8763.8     8829 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.6214  0.1441  0.3650  0.6075  1.2670 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  ID     (Intercept) 1.577    1.256   
## Number of obs: 8832, groups:  ID, 138
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.48452    0.11315  13.120  < 2e-16 ***
## cond_cc1    -0.20840    0.05335  -3.906 9.37e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## cond_cc1 -0.013

Robustness check 2: Learning models with the time covariate removed

Experiments 1 and 2

summary(m1_train_no_cov_regr <- glmer(correct ~  
                                      cond_cc + 
                                      (cond_cc |ID), 
          data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc + (cond_cc | ID)
##    Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  46182.9  46225.9 -23086.5  46172.9    39995 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3753 -1.1134  0.5036  0.6746  1.0940 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr 
##  ID     (Intercept) 0.3342   0.5781        
##         cond_cc1    0.4616   0.6794   -0.20
## Number of obs: 40000, groups:  ID, 125
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.97023    0.05308  18.279  < 2e-16 ***
## cond_cc1    -0.17612    0.06533  -2.696  0.00702 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## cond_cc1 -0.189
summary(m2_train_no_cov_regr <- glmer(correct ~  
                                      cond_cc + 
                                      (cond_cc | ID), 
          data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc + (cond_cc | ID)
##    Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  49610.0  49653.5 -24800.0  49600.0    44155 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.9719 -1.0789  0.4729  0.6802  1.0720 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr 
##  ID     (Intercept) 0.4287   0.6548        
##         cond_cc1    0.5173   0.7192   -0.09
## Number of obs: 44160, groups:  ID, 138
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.0489     0.0570  18.401  < 2e-16 ***
## cond_cc1     -0.1697     0.0657  -2.583  0.00981 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## cond_cc1 -0.088
sjPlot::tab_model(m1_train_full_regr)
  correct
Predictors Odds Ratios CI p
(Intercept) 2.83 2.51 – 3.20 <0.001
cond cc1 0.83 0.73 – 0.95 0.007
trial within condition 1.45 1.36 – 1.54 <0.001
Random Effects
σ2 3.29
τ00 ID 0.45
τ11 ID.cond_cc1 0.49
τ11 ID.scale(trial_within_condition) 0.12
ρ01 -0.19
0.82
ICC 0.17
N ID 125
Observations 40000
Marginal R2 / Conditional R2 0.035 / 0.202
sjPlot::tab_model(m2_train_full_regr)
  correct
Predictors Odds Ratios CI p
(Intercept) 3.16 2.75 – 3.63 <0.001
cond cc1 0.84 0.73 – 0.96 0.011
trial within condition 1.48 1.38 – 1.59 <0.001
Random Effects
σ2 3.29
τ00 ID 0.67
τ11 ID.cond_cc1 0.56
τ11 ID.scale(trial_within_condition) 0.17
ρ01 -0.07
0.87
ICC 0.23
N ID 138
Observations 44160
Marginal R2 / Conditional R2 0.037 / 0.257
sjPlot::tab_model(m1_test_full_regr)
  correct
Predictors Odds Ratios CI p
(Intercept) 4.80 3.80 – 6.07 <0.001
cond cc1 0.63 0.47 – 0.85 0.002
Random Effects
σ2 3.29
τ00 ID 1.56
τ11 ID.cond_cc1 1.97
ρ01 ID -0.17
ICC 0.38
N ID 125
Observations 8000
Marginal R2 / Conditional R2 0.010 / 0.390
sjPlot::tab_model(m2_test_full_regr)
  correct
Predictors Odds Ratios CI p
(Intercept) 5.32 4.18 – 6.78 <0.001
cond cc1 0.70 0.51 – 0.96 0.025
Random Effects
σ2 3.29
τ00 ID 1.85
τ11 ID.cond_cc1 2.47
ρ01 ID -0.15
ICC 0.43
N ID 138
Observations 8832
Marginal R2 / Conditional R2 0.005 / 0.432

Robustness check 3: Covariate for answer type

s1_train[s1_train$response=="sum", "resp_type_cc_cog"] <- .5
s1_train[s1_train$response=="difference", "resp_type_cc_cog"] <- -.5
s1_train[s1_train$response=="top", "resp_type_cc_cog"] <- 0
s1_train[s1_train$response=="bottom", "resp_type_cc_cog"] <- 0


s1_train[s1_train$response=="top", "resp_type_cc_overt"] <- .5
s1_train[s1_train$response=="bottom", "resp_type_cc_overt"] <- -.5
s1_train[s1_train$response=="sum", "resp_type_cc_overt"] <- 0
s1_train[s1_train$response=="difference", "resp_type_cc_overt"] <- 0

s1_sit[s1_sit$resp_as_category=="sum", "resp_type_cc_cog"] <- .5
s1_sit[s1_sit$resp_as_category=="difference", "resp_type_cc_cog"] <- -.5
s1_sit[s1_sit$resp_as_category=="top", "resp_type_cc_cog"] <- 0
s1_sit[s1_sit$resp_as_category=="bottom", "resp_type_cc_cog"] <- 0

s1_sit[s1_sit$resp_as_category=="top", "resp_type_cc_overt"] <- .5
s1_sit[s1_sit$resp_as_category=="bottom", "resp_type_cc_overt"] <- -.5
s1_sit[s1_sit$resp_as_category=="sum", "resp_type_cc_overt"] <- 0
s1_sit[s1_sit$resp_as_category=="difference", "resp_type_cc_overt"] <- 0

s2_train[s2_train$response=="alphabetize", "resp_type_cc_cog"] <- .5
s2_train[s2_train$response=="rev_alphabetize", "resp_type_cc_cog"] <- -.5
s2_train[s2_train$response=="slash", "resp_type_cc_cog"] <- 0
s2_train[s2_train$response=="backslash", "resp_type_cc_cog"] <- 0

s2_train[s2_train$response=="slash", "resp_type_cc_overt"] <- .5
s2_train[s2_train$response=="backslash", "resp_type_cc_overt"] <- -.5
s2_train[s2_train$response=="alphabetize", "resp_type_cc_overt"] <- 0
s2_train[s2_train$response=="rev_alphabetize", "resp_type_cc_overt"] <- 0

s2_sit[s2_sit$resp_as_category=="alphabetize", "resp_type_cc_cog"] <- .5
s2_sit[s2_sit$resp_as_category=="rev_alphabetize", "resp_type_cc_cog"] <- -.5
s2_sit[s2_sit$resp_as_category=="slash", "resp_type_cc_cog"] <- 0
s2_sit[s2_sit$resp_as_category=="backslash", "resp_type_cc_cog"] <- 0

s2_sit[s2_sit$resp_as_category=="slash", "resp_type_cc_overt"] <- .5
s2_sit[s2_sit$resp_as_category=="backslash", "resp_type_cc_overt"] <- -.5
s2_sit[s2_sit$resp_as_category=="alphabetize", "resp_type_cc_overt"] <- 0
s2_sit[s2_sit$resp_as_category=="rev_alphabetize", "resp_type_cc_overt"] <- 0

table(s1_train$resp_type_cc_cog)
## 
##  -0.5     0   0.5 
##  9020 20000 10980
table(s1_train$resp_type_cc_overt)
## 
##  -0.5     0   0.5 
##  9451 20000 10549
table(s1_sit$resp_type_cc_cog)
## 
## -0.5    0  0.5 
## 1851 4000 2149
table(s1_sit$resp_type_cc_overt)
## 
## -0.5    0  0.5 
## 1975 4000 2025
table(s2_train$resp_type_cc_cog)
## 
##  -0.5     0   0.5 
## 10182 22080 11898
table(s2_train$resp_type_cc_overt)
## 
##  -0.5     0   0.5 
## 11186 22080 10894
table(s2_sit$resp_type_cc_cog)
## 
## -0.5    0  0.5 
## 2024 4416 2392
table(s2_sit$resp_type_cc_overt)
## 
## -0.5    0  0.5 
## 2239 4416 2177

Separately coded contrast variables for cog and overt. Random slopes of these removed for test bc singular (which makes sense bc of few observations/subj)

summary(m1_train_full_regr_at <- glmer(correct ~  
                                      cond_cc + scale(trial_within_condition) + resp_type_cc_cog + resp_type_cc_overt +
                                      (cond_cc +  scale(trial_within_condition) + resp_type_cc_cog + resp_type_cc_overt |ID), 
          data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## correct ~ cond_cc + scale(trial_within_condition) + resp_type_cc_cog +  
##     resp_type_cc_overt + (cond_cc + scale(trial_within_condition) +  
##     resp_type_cc_cog + resp_type_cc_overt | ID)
##    Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    45040    45212   -22500    45000    39980 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -9.2192 -1.0643  0.4553  0.7031  1.4380 
## 
## Random effects:
##  Groups Name                          Variance Std.Dev. Corr                   
##  ID     (Intercept)                   0.44641  0.6681                          
##         cond_cc1                      0.49626  0.7045   -0.19                  
##         scale(trial_within_condition) 0.11598  0.3406    0.81  0.01            
##         resp_type_cc_cog              0.05453  0.2335    0.19 -0.31  0.33      
##         resp_type_cc_overt            0.07932  0.2816   -0.04  0.23  0.16  0.70
## Number of obs: 40000, groups:  ID, 125
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.06182    0.06142  17.287  < 2e-16 ***
## cond_cc1                      -0.18292    0.06840  -2.674  0.00749 ** 
## scale(trial_within_condition)  0.37205    0.03308  11.247  < 2e-16 ***
## resp_type_cc_cog              -0.07227    0.04151  -1.741  0.08165 .  
## resp_type_cc_overt            -0.08410    0.04532  -1.855  0.06353 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##              (Intr) cnd_c1 sc(__) rsp_typ_cc_c
## cond_cc1     -0.172                           
## scl(trl_w_)   0.753  0.004                    
## rsp_typ_cc_c  0.077 -0.175  0.164             
## rsp_typ_cc_v -0.034  0.137  0.088  0.197
# Singular 
# summary(m1_test_full_regr_at <- glmer(correct ~  cond_cc + resp_type_cc_cog + resp_type_cc_overt +
#                                         (cond_cc + resp_type_cc_cog + resp_type_cc_overt |ID), 
#           data=s1_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))

summary(m1_test_full_regr_at <- glmer(correct ~  cond_cc + resp_type_cc_cog + resp_type_cc_overt +
                                        (cond_cc  |ID), 
          data=s1_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## correct ~ cond_cc + resp_type_cc_cog + resp_type_cc_overt + (cond_cc |      ID)
##    Data: s1_sit
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   7857.3   7906.2  -3921.6   7843.3     7993 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0858  0.1296  0.2776  0.5823  1.2399 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr 
##  ID     (Intercept) 1.555    1.247         
##         cond_cc1    1.979    1.407    -0.17
## Number of obs: 8000, groups:  ID, 125
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         1.57061    0.11923  13.173  < 2e-16 ***
## cond_cc1           -0.45590    0.15119  -3.015  0.00257 ** 
## resp_type_cc_cog   -0.05969    0.08348  -0.715  0.47460    
## resp_type_cc_overt -0.07056    0.08935  -0.790  0.42972    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##              (Intr) cnd_c1 rsp_typ_cc_c
## cond_cc1     -0.175                    
## rsp_typ_cc_c -0.015 -0.021             
## rsp_typ_cc_v -0.007  0.011  0.000
summary(m2_train_full_regr_at <- glmer(correct ~  
                                      cond_cc + scale(trial_within_condition) + resp_type_cc_cog + resp_type_cc_overt +
                                      (cond_cc +  scale(trial_within_condition) + resp_type_cc_cog + resp_type_cc_overt |ID), 
          data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## correct ~ cond_cc + scale(trial_within_condition) + resp_type_cc_cog +  
##     resp_type_cc_overt + (cond_cc + scale(trial_within_condition) +  
##     resp_type_cc_cog + resp_type_cc_overt | ID)
##    Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  48297.6  48471.5 -24128.8  48257.6    44140 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -9.4618 -1.0313  0.4323  0.6859  1.4685 
## 
## Random effects:
##  Groups Name                          Variance Std.Dev. Corr                   
##  ID     (Intercept)                   0.6756   0.8220                          
##         cond_cc1                      0.5985   0.7736   -0.07                  
##         scale(trial_within_condition) 0.1666   0.4081    0.87 -0.04            
##         resp_type_cc_cog              0.2001   0.4474    0.01 -0.26  0.01      
##         resp_type_cc_overt            0.1730   0.4160   -0.27  0.00 -0.27  0.18
## Number of obs: 44160, groups:  ID, 138
## 
## Fixed effects:
##                                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.194005   0.071628  16.670   <2e-16 ***
## cond_cc1                      -0.163191   0.071098  -2.295   0.0217 *  
## scale(trial_within_condition)  0.398418   0.037258  10.693   <2e-16 ***
## resp_type_cc_cog              -0.122259   0.052879  -2.312   0.0208 *  
## resp_type_cc_overt            -0.001977   0.051645  -0.038   0.9695    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##              (Intr) cnd_c1 sc(__) rsp_typ_cc_c
## cond_cc1     -0.064                           
## scl(trl_w_)   0.821 -0.037                    
## rsp_typ_cc_c -0.001 -0.199  0.016             
## rsp_typ_cc_v -0.183  0.007 -0.170  0.087
# Singular  
# summary(m2_test_full_regr_at <- glmer(correct ~  cond_cc + resp_type_cc_cog + resp_type_cc_overt +
#                                         (cond_cc  + resp_type_cc_cog + resp_type_cc_overt  |ID), 
#           data=s2_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))

summary(m2_test_full_regr_at <- glmer(correct ~  cond_cc + resp_type_cc_cog + resp_type_cc_overt +
                                        (cond_cc  |ID), 
          data=s2_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## correct ~ cond_cc + resp_type_cc_cog + resp_type_cc_overt + (cond_cc |      ID)
##    Data: s2_sit
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   8427.9   8477.5  -4207.0   8413.9     8825 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.3845  0.1211  0.2303  0.5752  1.5798 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr 
##  ID     (Intercept) 1.848    1.359         
##         cond_cc1    2.469    1.571    -0.16
## Number of obs: 8832, groups:  ID, 138
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         1.67413    0.12343  13.564   <2e-16 ***
## cond_cc1           -0.35472    0.15863  -2.236   0.0253 *  
## resp_type_cc_cog   -0.05734    0.08376  -0.685   0.4936    
## resp_type_cc_overt  0.11865    0.08580   1.383   0.1667    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##              (Intr) cnd_c1 rsp_typ_cc_c
## cond_cc1     -0.159                    
## rsp_typ_cc_c -0.014 -0.022             
## rsp_typ_cc_v  0.010 -0.015  0.000

Robustness check 4: paired t-test on arcsine-transformed proportion data instead of mixed-effects model

s1_summs_id <- s1_train %>% group_by(condition, ID) %>% 
   summarize(m=mean(correct)) 
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
s2_summs_id <- s2_train %>% group_by(condition, ID) %>% 
   summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
s1_summs_test_id <- s1_sit %>% group_by(condition, ID) %>% 
   summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
s2_summs_test_id <- s2_sit %>% group_by(condition, ID) %>% 
   summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
t.test(asin(sqrt(m)) ~ condition, data=s1_summs_id, paired=TRUE)
## 
##  Paired t-test
## 
## data:  asin(sqrt(m)) by condition
## t = -2.5152, df = 124, p-value = 0.01318
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.061942537 -0.007385729
## sample estimates:
## mean difference 
##     -0.03466413
t.test(asin(sqrt(m)) ~ condition, data=s2_summs_id, paired=TRUE)
## 
##  Paired t-test
## 
## data:  asin(sqrt(m)) by condition
## t = -2.5139, df = 137, p-value = 0.0131
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.061281292 -0.007319414
## sample estimates:
## mean difference 
##     -0.03430035
t.test(asin(sqrt(m)) ~ condition, data=s1_summs_test_id, paired=TRUE)
## 
##  Paired t-test
## 
## data:  asin(sqrt(m)) by condition
## t = -2.6957, df = 124, p-value = 0.008001
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.12657517 -0.01939599
## sample estimates:
## mean difference 
##     -0.07298558
t.test(asin(sqrt(m)) ~ condition, data=s2_summs_test_id, paired=TRUE)
## 
##  Paired t-test
## 
## data:  asin(sqrt(m)) by condition
## t = -1.8641, df = 137, p-value = 0.06445
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.105733994  0.003119297
## sample estimates:
## mean difference 
##     -0.05130735

Additional check for resub: Is there moderation by difference in reaction times?

mean(s1_train[s1_train$condition=="cognitive", "rt_key_1"])-mean(s1_train[s1_train$condition=="overt", "rt_key_1"])
## [1] 246.5527
mean(s1_sit[s1_sit$condition=="cognitive", "rt_key_1"])-mean(s1_sit[s1_sit$condition=="overt", "rt_key_1"])
## [1] 188.8504
mean(s2_train[s2_train$condition=="cognitive", "rt_key_1"])-mean(s2_train[s2_train$condition=="overt", "rt_key_1"])
## [1] 680.3651
mean(s2_sit[s2_sit$condition=="cognitive", "rt_key_1"])-mean(s2_sit[s2_sit$condition=="overt", "rt_key_1"])
## [1] 643.3419
s1_tr_rt_inds <- s1_train %>% group_by(condition, ID) %>% summarize(m=mean(rt_key_1))   
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
s1_tr_diffs_inds <- 
  data.frame("diff"=s1_tr_rt_inds[s1_tr_rt_inds$condition=="cognitive", "m"]-s1_tr_rt_inds[s1_tr_rt_inds$condition=="overt", "m"],
             "ID"=s1_tr_rt_inds[s1_tr_rt_inds$condition=="cognitive", "ID"])
assert(all(s1_tr_rt_inds[s1_tr_rt_inds$condition=="cognitive", "ID"]==s1_tr_rt_inds[s1_tr_rt_inds$condition=="overt", "ID"]))

s1_te_rt_inds <- s1_sit %>% group_by(condition, ID) %>% summarize(m=mean(rt_key_1))   
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
s1_te_diffs_inds <- 
  data.frame("diff"=s1_te_rt_inds[s1_te_rt_inds$condition=="cognitive", "m"]-s1_te_rt_inds[s1_te_rt_inds$condition=="overt", "m"],
             "ID"=s1_te_rt_inds[s1_te_rt_inds$condition=="cognitive", "ID"])

s2_tr_rt_inds <- s2_train %>% group_by(condition, ID) %>% summarize(m=mean(rt_key_1))   
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
s2_tr_diffs_inds <- 
  data.frame("diff"=s2_tr_rt_inds[s2_tr_rt_inds$condition=="cognitive", "m"]-s2_tr_rt_inds[s2_tr_rt_inds$condition=="overt", "m"],
             "ID"=s2_tr_rt_inds[s2_tr_rt_inds$condition=="cognitive", "ID"])
assert(all(s2_tr_rt_inds[s2_tr_rt_inds$condition=="cognitive", "ID"]==s2_tr_rt_inds[s2_tr_rt_inds$condition=="overt", "ID"]))

s2_te_rt_inds <- s2_sit %>% group_by(condition, ID) %>% summarize(m=mean(rt_key_1))   
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
s2_te_diffs_inds <- 
  data.frame("diff"=s2_te_rt_inds[s2_te_rt_inds$condition=="cognitive", "m"]-s2_te_rt_inds[s2_te_rt_inds$condition=="overt", "m"],
             "ID"=s2_te_rt_inds[s2_te_rt_inds$condition=="cognitive", "ID"])

Add the RT diffs to dataframes

s1_IDs <- unique(s1_train$ID)
for (i in 1:length(s1_IDs)) {
  s1_train[s1_train$ID == s1_IDs[i], "rt_diff"] <- s1_tr_diffs_inds[s1_tr_diffs_inds$ID == s1_IDs[i], "m"]
  s1_sit[s1_sit$ID == s1_IDs[i], "rt_diff"] <- s1_te_diffs_inds[s1_te_diffs_inds$ID == s1_IDs[i], "m"]
}
# Spot checks 
# s1_train %>% filter(ID==25) %>% select(rt_diff)
# s1_tr_diffs_inds %>% filter(ID==25) %>% select(m)
# 
# s1_sit %>% filter(ID==25) %>% select(rt_diff)
# s1_te_diffs_inds %>% filter(ID==25) %>% select(m)
# 
# s1_train %>% filter(ID==102) %>% select(rt_diff)
# s1_tr_diffs_inds %>% filter(ID==102) %>% select(m)
# 
# s1_sit %>% filter(ID==102) %>% select(rt_diff)
# s1_te_diffs_inds %>% filter(ID==102) %>% select(m)
s2_IDs <- unique(s2_train$ID)
for (i in 1:length(s1_IDs)) {
  s2_train[s2_train$ID == s2_IDs[i], "rt_diff"] <- s2_tr_diffs_inds[s2_tr_diffs_inds$ID == s2_IDs[i], "m"]
  s2_sit[s2_sit$ID == s2_IDs[i], "rt_diff"] <- s2_te_diffs_inds[s2_te_diffs_inds$ID == s2_IDs[i], "m"]
}
# Spot checks 
# s2_train %>% filter(ID==25) %>% select(rt_diff)
# s2_tr_diffs_inds %>% filter(ID==25) %>% select(m)
# 
# s2_sit %>% filter(ID==25) %>% select(rt_diff)
# s2_te_diffs_inds %>% filter(ID==25) %>% select(m)
# 
# s2_train %>% filter(ID==102) %>% select(rt_diff)
# s2_tr_diffs_inds %>% filter(ID==102) %>% select(m)
# # 
# s2_sit %>% filter(ID==102) %>% select(rt_diff)
# s2_te_diffs_inds %>% filter(ID==102) %>% select(m)

Examine as potential moderator

Reduced complexity where needed because more complex models were singular

summary(m1_train_full_regr <- glmer(correct ~  
                                      cond_cc*scale(rt_diff) + scale(trial_within_condition) + 
                                      (cond_cc + scale(trial_within_condition)  |ID), 
          data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc * scale(rt_diff) + scale(trial_within_condition) +  
##     (cond_cc + scale(trial_within_condition) | ID)
##    Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  45057.4  45152.0 -22517.7  45035.4    39989 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -8.8320 -1.0683  0.4574  0.7050  1.2730 
## 
## Random effects:
##  Groups Name                          Variance Std.Dev. Corr       
##  ID     (Intercept)                   0.4451   0.6672              
##         cond_cc1                      0.4900   0.7000   -0.18      
##         scale(trial_within_condition) 0.1156   0.3400    0.82  0.01
## Number of obs: 40000, groups:  ID, 125
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.04184    0.06112  17.046  < 2e-16 ***
## cond_cc1                      -0.18034    0.06728  -2.681  0.00735 ** 
## scale(rt_diff)                 0.02641    0.03960   0.667  0.50490    
## scale(trial_within_condition)  0.36992    0.03300  11.209  < 2e-16 ***
## cond_cc1:scale(rt_diff)       -0.04805    0.06688  -0.718  0.47249    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cnd_c1 scl(_) sc(__)
## cond_cc1    -0.170                     
## scl(rt_dff)  0.000  0.000              
## scl(trl_w_)  0.759  0.010  0.000       
## cnd_cc1:(_)  0.000 -0.001 -0.274  0.000
summary(m1_test_full_regr <- glmer(correct ~  cond_cc*scale(rt_diff) + (cond_cc |ID), 
          data=s1_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc * scale(rt_diff) + (cond_cc | ID)
##    Data: s1_sit
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   7856.2   7905.1  -3921.1   7842.2     7993 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0633  0.1280  0.2738  0.5815  1.2273 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr 
##  ID     (Intercept) 1.537    1.240         
##         cond_cc1    1.945    1.395    -0.15
## Number of obs: 8000, groups:  ID, 125
## 
## Fixed effects:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               1.5680     0.1186  13.221  < 2e-16 ***
## cond_cc1                 -0.4534     0.1502  -3.018  0.00254 ** 
## scale(rt_diff)            0.1212     0.1159   1.046  0.29555    
## cond_cc1:scale(rt_diff)  -0.1729     0.1412  -1.225  0.22051    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cnd_c1 scl(_)
## cond_cc1    -0.160              
## scl(rt_dff)  0.006 -0.015       
## cnd_cc1:(_) -0.016  0.015 -0.145

Reduced complexity for these two because more complex models singular

summary(m2_train_full_regr <- glmer(correct ~  
                                      cond_cc*scale(rt_diff) + scale(trial_within_condition) + 
                                      (cond_cc + scale(trial_within_condition)  |ID), 
          data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc * scale(rt_diff) + scale(trial_within_condition) +  
##     (cond_cc + scale(trial_within_condition) | ID)
##    Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  43489.9  43584.5 -21734.0  43467.9    39989 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -10.2304  -1.0286   0.4327   0.6848   1.4586 
## 
## Random effects:
##  Groups Name                          Variance Std.Dev. Corr       
##  ID     (Intercept)                   0.7250   0.8515              
##         cond_cc1                      0.5118   0.7154   -0.08      
##         scale(trial_within_condition) 0.1803   0.4247    0.88 -0.07
## Number of obs: 40000, groups:  ID, 125
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.17379    0.07769  15.109  < 2e-16 ***
## cond_cc1                      -0.15927    0.06908  -2.306  0.02113 *  
## scale(rt_diff)                 0.05355    0.04297   1.246  0.21271    
## scale(trial_within_condition)  0.40406    0.04059   9.954  < 2e-16 ***
## cond_cc1:scale(rt_diff)       -0.22740    0.06887  -3.302  0.00096 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cnd_c1 scl(_) sc(__)
## cond_cc1    -0.074                     
## scl(rt_dff)  0.001 -0.007              
## scl(trl_w_)  0.832 -0.064 -0.001       
## cnd_cc1:(_) -0.006  0.006 -0.032 -0.002
summary(m2_test_full_regr <- glmer(correct ~  cond_cc*scale(rt_diff) + (cond_cc |ID), 
          data=s2_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc * scale(rt_diff) + (cond_cc | ID)
##    Data: s2_sit
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   7490.3   7539.2  -3738.2   7476.3     7993 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.6644  0.1168  0.2216  0.5651  1.5704 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr 
##  ID     (Intercept) 1.904    1.380         
##         cond_cc1    2.307    1.519    -0.09
## Number of obs: 8000, groups:  ID, 125
## 
## Fixed effects:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               1.7299     0.1318  13.124   <2e-16 ***
## cond_cc1                 -0.2647     0.1640  -1.614   0.1066    
## scale(rt_diff)            0.2416     0.1304   1.852   0.0640 .  
## cond_cc1:scale(rt_diff)  -0.3419     0.1580  -2.164   0.0305 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cnd_c1 scl(_)
## cond_cc1    -0.093              
## scl(rt_dff)  0.020 -0.018       
## cnd_cc1:(_) -0.017  0.052 -0.100

Additional check for resub 2: Is performance better when the first learning block of the condition happens second, reflective of meta-learning

s1_with_cond_order <- read.csv("../data/cleaned_files/s1_train_deident_with_first_cond.csv")
s2_with_cond_order <- read.csv("../data/cleaned_files/s2_train_with_exclusions_IDed_with_first_cond_wdeID.csv")

Put the first assignments into regular dfs

unique_IDs <- unique(s1_train$ID)

for (s in 1:length(unique_IDs)) {
  s1_train[s1_train$ID == unique_IDs[s], "first_cond"] <- 
    s1_with_cond_order[s1_with_cond_order$deident_ID == unique_IDs[s], "first_cond"]
  s1_sit[s1_sit$ID == unique_IDs[s], "first_cond"] <- 
    unique(s1_with_cond_order[s1_with_cond_order$deident_ID == unique_IDs[s], "first_cond"])
}
s1_train %>% group_by(condition, first_cond) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## # A tibble: 4 × 3
## # Groups:   condition [2]
##   condition first_cond        m
##   <chr>     <chr>         <dbl>
## 1 cognitive cognTraining  0.670
## 2 cognitive overtTraining 0.715
## 3 overt     cognTraining  0.723
## 4 overt     overtTraining 0.716
s1_sit %>% group_by(condition, first_cond) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## # A tibble: 4 × 3
## # Groups:   condition [2]
##   condition first_cond        m
##   <chr>     <chr>         <dbl>
## 1 cognitive cognTraining  0.703
## 2 cognitive overtTraining 0.752
## 3 overt     cognTraining  0.790
## 4 overt     overtTraining 0.765
unique_IDs <- unique(s2_train$ID)

for (s in 1:length(unique_IDs)) {
  s2_train[s2_train$ID == unique_IDs[s], "first_cond"] <- 
    s2_with_cond_order[s2_with_cond_order$de_ID == unique_IDs[s], "first_cond"]
  s2_sit[s2_sit$ID == unique_IDs[s], "first_cond"] <- 
    unique(s2_with_cond_order[s2_with_cond_order$de_ID == unique_IDs[s], "first_cond"])
}
s2_train %>% group_by(condition, first_cond) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## # A tibble: 4 × 3
## # Groups:   condition [2]
##   condition first_cond     m
##   <chr>     <chr>      <dbl>
## 1 cognitive cognitive  0.667
## 2 cognitive overt      0.733
## 3 overt     cognitive  0.734
## 4 overt     overt      0.728
s2_sit %>% group_by(condition, first_cond) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## # A tibble: 4 × 3
## # Groups:   condition [2]
##   condition first_cond     m
##   <chr>     <chr>      <dbl>
## 1 cognitive cognitive  0.687
## 2 cognitive overt      0.784
## 3 overt     cognitive  0.775
## 4 overt     overt      0.769

Random assignment of first condition was largely even, with more assignments to cognitive first in study 1 and overt first in condition 2

table(s1_train$first_cond)
## 
##  cognTraining overtTraining 
##         20480         19520
table(s2_train$first_cond)
## 
## cognitive     overt 
##     20480     23680
# Specify overt as baseline factor so that negative effect corresponds to worst performance  
s1_train$fc_cond <- factor(s1_train$first_cond, levels=c("overtTraining", "cognTraining"))
contrasts(s1_train$fc_cond) <- c(-.5, .5)
#head(s1_train$fc_cond)

s2_train$fc_cond <- factor(s2_train$first_cond, levels=c("overt", "cognitive"))
contrasts(s2_train$fc_cond) <- c(-.5, .5)

s1_sit$fc_cond <- factor(s1_sit$first_cond, levels=c("overtTraining", "cognTraining"))
contrasts(s1_sit$fc_cond) <- c(-.5, .5)

s2_sit$fc_cond <- factor(s2_sit$first_cond, levels=c("overt", "cognitive"))
contrasts(s2_sit$fc_cond) <- c(-.5, .5)

Consistent pattern whereby the condition difference is most pronounced when cognitive is first; main effect remains significant

# Reduced RE from interaction bc of non-convergence with full  
summary(m1_train_full_regr <- glmer(correct ~  
                                      cond_cc*fc_cond + scale(trial_within_condition) + 
                                      (cond_cc + fc_cond +  scale(trial_within_condition) |ID), 
          data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc * fc_cond + scale(trial_within_condition) +  
##     (cond_cc + fc_cond + scale(trial_within_condition) | ID)
##    Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  45053.5  45182.5 -22511.8  45023.5    39985 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -8.9434 -1.0673  0.4574  0.7052  1.2873 
## 
## Random effects:
##  Groups Name                          Variance Std.Dev. Corr             
##  ID     (Intercept)                   0.3225   0.5679                    
##         cond_cc1                      0.4698   0.6854   -0.27            
##         fc_cond1                      0.5118   0.7154   -0.30  0.00      
##         scale(trial_within_condition) 0.1154   0.3398    0.96 -0.02 -0.22
## Number of obs: 40000, groups:  ID, 125
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.03334    0.06137  16.839  < 2e-16 ***
## cond_cc1                      -0.17725    0.06609  -2.682  0.00732 ** 
## fc_cond1                      -0.02511    0.08103  -0.310  0.75662    
## scale(trial_within_condition)  0.36981    0.03298  11.212  < 2e-16 ***
## cond_cc1:fc_cond1             -0.29105    0.13287  -2.191  0.02849 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cnd_c1 fc_cn1 sc(__)
## cond_cc1    -0.210                     
## fc_cond1    -0.285  0.006              
## scl(trl_w_)  0.758 -0.015 -0.152       
## cnd_cc1:f_1  0.001 -0.025 -0.301 -0.003
# Reduced for same reason  
summary(m1_test_full_regr <- glmer(correct ~  cond_cc*fc_cond + (cond_cc + fc_cond|ID), 
          data=s1_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc * fc_cond + (cond_cc + fc_cond | ID)
##    Data: s1_sit
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   7847.1   7917.0  -3913.6   7827.1     7990 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.2749  0.1146  0.2742  0.5759  1.2486 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr       
##  ID     (Intercept) 1.197    1.094               
##         cond_cc1    2.004    1.416    -0.14      
##         fc_cond1    1.889    1.374    -0.23 -0.64
## Number of obs: 8000, groups:  ID, 125
## 
## Fixed effects:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         1.5963     0.1240  12.874  < 2e-16 ***
## cond_cc1           -0.4114     0.1520  -2.707  0.00679 ** 
## fc_cond1           -0.1968     0.2456  -0.801  0.42295    
## cond_cc1:fc_cond1  -0.9458     0.3070  -3.081  0.00207 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cnd_c1 fc_cn1
## cond_cc1    -0.113              
## fc_cond1    -0.224 -0.299       
## cnd_cc1:f_1 -0.314 -0.063 -0.106
# Run with full interaction because 
summary(m2_train_full_regr <- glmer(correct ~  
                                      cond_cc*fc_cond + scale(trial_within_condition) + 
                                      (cond_cc*fc_cond +  scale(trial_within_condition) | ID), 
          data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc * fc_cond + scale(trial_within_condition) +  
##     (cond_cc * fc_cond + scale(trial_within_condition) | ID)
##    Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  48362.6  48536.6 -24161.3  48322.6    44140 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -10.2482  -1.0364   0.4408   0.6878   1.4992 
## 
## Random effects:
##  Groups Name                          Variance Std.Dev. Corr                   
##  ID     (Intercept)                   0.5537   0.7441                          
##         cond_cc1                      0.3106   0.5573   -0.28                  
##         fc_cond1                      0.4216   0.6493   -0.11 -0.30            
##         scale(trial_within_condition) 0.1689   0.4110    0.94 -0.19 -0.16      
##         cond_cc1:fc_cond1             0.8624   0.9287   -0.29 -0.10  0.14 -0.23
## Number of obs: 44160, groups:  ID, 138
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.14228    0.07120  16.042  < 2e-16 ***
## cond_cc1                      -0.21284    0.06739  -3.158  0.00159 ** 
## fc_cond1                       0.02138    0.08324   0.257  0.79734    
## scale(trial_within_condition)  0.39290    0.03745  10.491  < 2e-16 ***
## cond_cc1:fc_cond1             -0.52759    0.13473  -3.916 9.01e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cnd_c1 fc_cn1 sc(__)
## cond_cc1    -0.179                     
## fc_cond1     0.007 -0.254              
## scl(trl_w_)  0.813 -0.133 -0.088       
## cnd_cc1:f_1 -0.262 -0.001 -0.109 -0.132
summary(m2_test_full_regr <- glmer(correct ~  cond_cc*fc_cond + (cond_cc*fc_cond|ID), 
          data=s2_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ cond_cc * fc_cond + (cond_cc * fc_cond | ID)
##    Data: s2_sit
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   8430.5   8529.7  -4201.3   8402.5     8818 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.2679  0.1224  0.2404  0.5711  1.5615 
## 
## Random effects:
##  Groups Name              Variance Std.Dev. Corr             
##  ID     (Intercept)       1.411    1.188                     
##         cond_cc1          1.957    1.399    -0.21            
##         fc_cond1          1.651    1.285    -0.28 -0.07      
##         cond_cc1:fc_cond1 1.458    1.208    -0.13 -0.05 -0.10
## Number of obs: 8832, groups:  ID, 138
## 
## Fixed effects:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         1.6640     0.1221  13.625  < 2e-16 ***
## cond_cc1           -0.3674     0.1554  -2.364  0.01808 *  
## fc_cond1           -0.4580     0.2441  -1.876  0.06068 .  
## cond_cc1:fc_cond1  -0.8071     0.3107  -2.598  0.00939 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cnd_c1 fc_cn1
## cond_cc1    -0.183              
## fc_cond1    -0.158 -0.098       
## cnd_cc1:f_1 -0.098 -0.007 -0.183

For resub, correlation between training and test performance

cor.test(
  as.numeric(unlist(s1_train_summs[s1_train_summs$condition=="overt", "m"])),
  as.numeric(unlist(s1_test_summs[s1_test_summs$condition=="overt", "m"]))
)
## 
##  Pearson's product-moment correlation
## 
## data:  as.numeric(unlist(s1_train_summs[s1_train_summs$condition == "overt", "m"])) and as.numeric(unlist(s1_test_summs[s1_test_summs$condition == "overt", "m"]))
## t = 12.528, df = 123, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6599115 0.8169410
## sample estimates:
##       cor 
## 0.7487498
cor.test(
  as.numeric(unlist(s1_train_summs[s1_train_summs$condition=="cognitive", "m"])),
  as.numeric(unlist(s1_test_summs[s1_test_summs$condition=="cognitive", "m"]))
)
## 
##  Pearson's product-moment correlation
## 
## data:  as.numeric(unlist(s1_train_summs[s1_train_summs$condition == "cognitive", "m"])) and as.numeric(unlist(s1_test_summs[s1_test_summs$condition == "cognitive", "m"]))
## t = 14.454, df = 123, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7177575 0.8504835
## sample estimates:
##       cor 
## 0.7933662
cor.test(
  as.numeric(unlist(s2_train_summs[s2_train_summs$condition=="overt", "m"])),
  as.numeric(unlist(s2_test_summs[s2_test_summs$condition=="overt", "m"]))
)
## 
##  Pearson's product-moment correlation
## 
## data:  as.numeric(unlist(s2_train_summs[s2_train_summs$condition == "overt", "m"])) and as.numeric(unlist(s2_test_summs[s2_test_summs$condition == "overt", "m"]))
## t = 11.385, df = 136, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6016955 0.7751783
## sample estimates:
##       cor 
## 0.6985625
cor.test(
  as.numeric(unlist(s2_train_summs[s2_train_summs$condition=="cognitive", "m"])),
  as.numeric(unlist(s2_test_summs[s2_test_summs$condition=="cognitive", "m"]))
)
## 
##  Pearson's product-moment correlation
## 
## data:  as.numeric(unlist(s2_train_summs[s2_train_summs$condition == "cognitive", "m"])) and as.numeric(unlist(s2_test_summs[s2_test_summs$condition == "cognitive", "m"]))
## t = 13.631, df = 136, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6789692 0.8225245
## sample estimates:
##       cor 
## 0.7598609

Delay effects

Create sum contrast codes for delay

# Specify overt as baseline factor so that negative effect corresponds to worst performance  
s1_train$cond_cc <- factor(s1_train$condition, levels=c("overt", "cognitive"))
contrasts(s1_train$cond_cc) <- c(-.5, .5)

s2_train$cond_cc <- factor(s2_train$condition, levels=c("overt", "cognitive"))
contrasts(s2_train$cond_cc) <- c(-.5, .5)
#head(factor(if_else(s1_train$delay==0, "no_delay", "delay"), levels=c("no_delay", "delay")))
s1_train$delay_cc <- factor(if_else(s1_train$delay==0, "no_delay", "delay"), levels=c("no_delay", "delay"))
contrasts(s1_train$delay_cc) <- c(-.5, .5)

s2_train$delay_cc <- factor(if_else(s2_train$delay==0, "no_delay", "delay"), levels=c("no_delay", "delay"))
contrasts(s2_train$delay_cc) <- c(-.5, .5)

Building up in complexity…

Effect of delay with no moderation

summary(m1_delay <- glmer(correct ~  delay_cc + (delay_cc|ID), 
          data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ delay_cc + (delay_cc | ID)
##    Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  43647.9  43690.6 -21818.9  43637.9    37995 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2305 -1.1242  0.5124  0.6818  1.0042 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr
##  ID     (Intercept) 0.32784  0.5726       
##         delay_cc1   0.02447  0.1564   0.43
## Number of obs: 38000, groups:  ID, 125
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.07222    0.05324   20.14   <2e-16 ***
## delay_cc1   -0.34260    0.03240  -10.57   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## delay_cc1 0.043
summary(m2_delay <- glmer(correct ~  delay_cc + (delay_cc|ID), 
          data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ delay_cc + (delay_cc | ID)
##    Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  47084.2  47127.5 -23537.1  47074.2    41947 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.2244 -1.1013  0.4908  0.6573  1.0211 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr
##  ID     (Intercept) 0.42730  0.6537       
##         delay_cc1   0.01198  0.1095   0.44
## Number of obs: 41952, groups:  ID, 138
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.13294    0.05749   19.71   <2e-16 ***
## delay_cc1   -0.30347    0.03024  -10.04   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## delay_cc1 -0.001

Now moderated by condition

#Singular for both 
# summary(m1_delay_cond <- glmer(correct ~  delay_cc*cond_cc + (delay_cc*cond_cc|ID), 
#           data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
# 
# summary(m2_delay_cond <- glmer(correct ~  delay_cc*cond_cc + (delay_cc*cond_cc|ID), 
#           data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))

# Simplified  
summary(m1_delay_cond <- glmer(correct ~  delay_cc*cond_cc + (delay_cc + cond_cc|ID),
          data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ delay_cc * cond_cc + (delay_cc + cond_cc | ID)
##    Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  43090.3  43175.7 -21535.1  43070.3    37990 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0844 -1.0664  0.4738  0.6770  1.1528 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr       
##  ID     (Intercept) 0.35727  0.5977              
##         delay_cc1   0.02429  0.1559    0.37      
##         cond_cc1    0.50501  0.7106   -0.22 -0.17
## Number of obs: 38000, groups:  ID, 125
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         1.10896    0.05554  19.967   <2e-16 ***
## delay_cc1          -0.35536    0.03271 -10.864   <2e-16 ***
## cond_cc1           -0.17949    0.07027  -2.554   0.0106 *  
## delay_cc1:cond_cc1 -0.01904    0.05747  -0.331   0.7404    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) dly_c1 cnd_c1
## delay_cc1    0.019              
## cond_cc1    -0.196 -0.055       
## dly_cc1:c_1  0.006 -0.041 -0.226
summary(m2_delay_cond <- glmer(correct ~  delay_cc*cond_cc + (delay_cc + cond_cc|ID),
          data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ delay_cc * cond_cc + (delay_cc + cond_cc | ID)
##    Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  46388.5  46474.9 -23184.2  46368.5    41942 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.9099 -1.0345  0.4587  0.6672  1.1512 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr       
##  ID     (Intercept) 0.45639  0.6756              
##         delay_cc1   0.01049  0.1024    0.55      
##         cond_cc1    0.56015  0.7484   -0.09 -0.21
## Number of obs: 41952, groups:  ID, 138
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         1.17084    0.05936  19.725   <2e-16 ***
## delay_cc1          -0.30464    0.03032 -10.046   <2e-16 ***
## cond_cc1           -0.13973    0.07005  -1.995   0.0461 *  
## delay_cc1:cond_cc1 -0.13157    0.05525  -2.381   0.0172 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) dly_c1 cnd_c1
## delay_cc1    0.023              
## cond_cc1    -0.085 -0.054       
## dly_cc1:c_1  0.001 -0.025 -0.214

Model including appropriate covariates at the highest complexity that would fit (removed trial-within-cond random slope and delay*cond slope due to singular errors in studies 2 and 1 respectively)

# summary(m1_full_delay_cond <- glmer(correct ~  delay_cc*cond_cc + scale(trial_within_condition) +
#                                  (delay_cc + cond_cc + scale(trial_within_condition)|ID),
#           data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
# 
# # Singular 
# summary(m2_full_delay_cond <- glmer(correct ~  delay_cc*cond_cc + scale(trial_within_condition) +
#                                  (delay_cc + cond_cc + scale(trial_within_condition)|ID),
#           data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))

# Also singular  
# summary(m1_full_delay_cond <- glmer(correct ~  delay_cc*cond_cc + scale(trial_within_condition) +
#                                  (delay_cc*cond_cc + 1|ID),
#           data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))

summary(m1_full_delay_cond <- glmer(correct ~  delay_cc*cond_cc + scale(trial_within_condition) +
                                 (delay_cc + cond_cc |ID),
          data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ delay_cc * cond_cc + scale(trial_within_condition) +  
##     (delay_cc + cond_cc | ID)
##    Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  42576.5  42670.5 -21277.2  42554.5    37989 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.2409 -1.0187  0.4627  0.6726  1.4465 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr       
##  ID     (Intercept) 0.36544  0.6045              
##         delay_cc1   0.02941  0.1715    0.36      
##         cond_cc1    0.51948  0.7208   -0.21 -0.18
## Number of obs: 38000, groups:  ID, 125
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.11508    0.05619  19.845  < 2e-16 ***
## delay_cc1                     -0.36388    0.03352 -10.857  < 2e-16 ***
## cond_cc1                      -0.18364    0.07121  -2.579  0.00991 ** 
## scale(trial_within_condition)  0.27402    0.01215  22.553  < 2e-16 ***
## delay_cc1:cond_cc1            -0.01022    0.05790  -0.177  0.85983    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) dly_c1 cnd_c1 sc(__)
## delay_cc1    0.029                     
## cond_cc1    -0.194 -0.064              
## scl(trl_w_)  0.018 -0.019 -0.004       
## dly_cc1:c_1  0.007 -0.041 -0.225  0.006
summary(m2_full_delay_cond <- glmer(correct ~  delay_cc*cond_cc + scale(trial_within_condition) +
                                 (delay_cc + cond_cc|ID),
          data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ delay_cc * cond_cc + scale(trial_within_condition) +  
##     (delay_cc + cond_cc | ID)
##    Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  45878.5  45973.6 -22928.2  45856.5    41941 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.8197 -0.9779  0.4437  0.6526  1.4299 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr       
##  ID     (Intercept) 0.46804  0.6841              
##         delay_cc1   0.01325  0.1151    0.47      
##         cond_cc1    0.57474  0.7581   -0.09 -0.18
## Number of obs: 41952, groups:  ID, 138
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.17603    0.06016  19.549   <2e-16 ***
## delay_cc1                     -0.30940    0.03083 -10.036   <2e-16 ***
## cond_cc1                      -0.14542    0.07100  -2.048   0.0405 *  
## scale(trial_within_condition)  0.26368    0.01172  22.489   <2e-16 ***
## delay_cc1:cond_cc1            -0.11654    0.05571  -2.092   0.0364 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) dly_c1 cnd_c1 sc(__)
## delay_cc1    0.016                     
## cond_cc1    -0.084 -0.052              
## scl(trl_w_)  0.016 -0.013 -0.004       
## dly_cc1:c_1  0.001 -0.025 -0.214  0.010
car::vif(m1_full_delay_cond)
##                      delay_cc                       cond_cc 
##                      1.007763                      1.059556 
## scale(trial_within_condition)              delay_cc:cond_cc 
##                      1.000389                      1.056908
car::vif(m2_full_delay_cond)
##                      delay_cc                       cond_cc 
##                      1.004255                      1.051818 
## scale(trial_within_condition)              delay_cc:cond_cc 
##                      1.000268                      1.049647

For resub, rerun interaction with continuous delay

summary(m1_full_delay_cond_cont <- glmer(correct ~  scale(delay)*cond_cc + scale(trial_within_condition) +
                                 (cond_cc|ID),
          data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ scale(delay) * cond_cc + scale(trial_within_condition) +  
##     (cond_cc | ID)
##    Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  42668.7  42737.1 -21326.4  42652.7    37992 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4888 -1.0349  0.4694  0.6753  1.5549 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr 
##  ID     (Intercept) 0.3873   0.6224        
##         cond_cc1    0.5156   0.7181   -0.22
## Number of obs: 38000, groups:  ID, 125
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.02230    0.05713  17.894  < 2e-16 ***
## scale(delay)                  -0.11266    0.01154  -9.766  < 2e-16 ***
## cond_cc1                      -0.18989    0.06914  -2.746  0.00603 ** 
## scale(trial_within_condition)  0.27169    0.01212  22.421  < 2e-16 ***
## scale(delay):cond_cc1          0.02219    0.02307   0.962  0.33610    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) scl(d) cnd_c1 sc(__)
## scale(dely) -0.009                     
## cond_cc1    -0.208  0.003              
## scl(trl_w_)  0.015 -0.015 -0.003       
## scl(dly):_1  0.002 -0.018 -0.015 -0.001
# Singular 
# summary(m1_full_delay_cond_cont <- glmer(correct ~  scale(delay)*cond_cc + scale(trial_within_condition) +
#                                  (scale(delay) + cond_cc|ID),
#           data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
# Singular 
# summary(m1_full_delay_cond_cont <- glmer(correct ~  scale(delay)*cond_cc + scale(trial_within_condition) +
#                                  (scale(delay)*cond_cc |ID),
#           data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
summary(m2_full_delay_cond_cont <- glmer(correct ~  scale(delay)*cond_cc + scale(trial_within_condition) +
                                 ( cond_cc|ID),
          data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: correct ~ scale(delay) * cond_cc + scale(trial_within_condition) +  
##     (cond_cc | ID)
##    Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  45932.2  46001.3 -22958.1  45916.2    41944 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.8928 -0.9866  0.4450  0.6548  1.5689 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr 
##  ID     (Intercept) 0.4873   0.6981        
##         cond_cc1    0.5749   0.7582   -0.10
## Number of obs: 41952, groups:  ID, 138
## 
## Fixed effects:
##                                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.097613   0.060775  18.060   <2e-16 ***
## scale(delay)                  -0.108518   0.011140  -9.741   <2e-16 ***
## cond_cc1                      -0.177724   0.069359  -2.562   0.0104 *  
## scale(trial_within_condition)  0.263420   0.011711  22.494   <2e-16 ***
## scale(delay):cond_cc1         -0.004329   0.022280  -0.194   0.8459    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) scl(d) cnd_c1 sc(__)
## scale(dely) -0.009                     
## cond_cc1    -0.092  0.001              
## scl(trl_w_)  0.014 -0.010 -0.003       
## scl(dly):_1  0.001 -0.030 -0.015 -0.003
# Also singular
# summary(m2_full_delay_cond_cont <- glmer(correct ~  scale(delay)*cond_cc + scale(trial_within_condition) +
#                                  (scale(delay) + cond_cc|ID),
#           data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))

# # singular
# summary(m2_full_delay_cond_cont <- glmer(correct ~  scale(delay)*cond_cc + scale(trial_within_condition) +
#                                  (scale(delay)*cond_cc|ID),
#           data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))

Heterogeneity in condition effect

Calculate some bxal differences

Study 1

s1_perf <- s1_train %>% group_by(condition, ID) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
o_min_c_perf_s1 <- as.numeric(unlist(s1_perf[s1_perf$condition=="overt","m"]-
                    s1_perf[s1_perf$condition=="cognitive", "m"]))

# .. and test subject-level differences in performance  
s1_perf_test <- s1_sit %>% group_by(condition, ID) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
o_min_c_testperf_s1 <- as.numeric(unlist(s1_perf_test[s1_perf_test$condition=="overt","m"]-
                    s1_perf_test[s1_perf_test$condition=="cognitive", "m"]))

s1_perf_overall <- s1_train %>% group_by(ID) %>% summarize(m=mean(correct))
s1_test_perf_overall <- s1_sit %>% group_by(ID) %>% summarize(m=mean(correct))

# Modest correlations across conditions 
cor.test(unlist(s1_perf[s1_perf$condition=="overt","m"]), unlist(s1_perf[s1_perf$condition=="cognitive","m"]))
## 
##  Pearson's product-moment correlation
## 
## data:  unlist(s1_perf[s1_perf$condition == "overt", "m"]) and unlist(s1_perf[s1_perf$condition == "cognitive", "m"])
## t = 5.461, df = 123, p-value = 2.507e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2885290 0.5729171
## sample estimates:
##       cor 
## 0.4417538
cor.test(unlist(s1_perf_test[s1_perf_test$condition=="overt","m"]), unlist(s1_perf_test[s1_perf_test$condition=="cognitive","m"]))
## 
##  Pearson's product-moment correlation
## 
## data:  unlist(s1_perf_test[s1_perf_test$condition == "overt", "m"]) and unlist(s1_perf_test[s1_perf_test$condition == "cognitive", "m"])
## t = 5.6112, df = 123, p-value = 1.263e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2995928 0.5809966
## sample estimates:
##       cor 
## 0.4514492

Examining if between-condition diffs are more common at the highest level of perf, because at lower ends more noise in bx. But not much relationship

cor.test(s1_perf_overall$m, o_min_c_perf_s1)
## 
##  Pearson's product-moment correlation
## 
## data:  s1_perf_overall$m and o_min_c_perf_s1
## t = 1.0074, df = 123, p-value = 0.3157
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.08652362  0.26190545
## sample estimates:
##        cor 
## 0.09045834
ComparePars(s1_perf_overall$m, o_min_c_perf_s1, use_identity_line = 0)

Study 2

s2_perf <- s2_train %>% group_by(condition, ID) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
o_min_c_perf_s2 <- as.numeric(unlist(s2_perf[s2_perf$condition=="overt","m"]-
                    s2_perf[s2_perf$condition=="cognitive", "m"]))

s2_perf_test <- s2_sit %>% group_by(condition, ID) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
o_min_c_testperf_s2 <- as.numeric(unlist(s2_perf_test[s2_perf_test$condition=="overt","m"]-
                    s2_perf_test[s2_perf_test$condition=="cognitive", "m"]))


s2_perf_overall <- s2_train %>% group_by(ID) %>% summarize(m=mean(correct))
s2_test_perf_overall <- s2_sit %>% group_by(ID) %>% summarize(m=mean(correct))

cor.test(unlist(s2_perf[s2_perf$condition=="overt","m"]), unlist(s2_perf[s2_perf$condition=="cognitive","m"]))
## 
##  Pearson's product-moment correlation
## 
## data:  unlist(s2_perf[s2_perf$condition == "overt", "m"]) and unlist(s2_perf[s2_perf$condition == "cognitive", "m"])
## t = 6.0926, df = 136, p-value = 1.074e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3207662 0.5848974
## sample estimates:
##       cor 
## 0.4630508
cor.test(unlist(s2_perf_test[s2_perf_test$condition=="overt","m"]), unlist(s2_perf_test[s2_perf_test$condition=="cognitive","m"]))
## 
##  Pearson's product-moment correlation
## 
## data:  unlist(s2_perf_test[s2_perf_test$condition == "overt", "m"]) and unlist(s2_perf_test[s2_perf_test$condition == "cognitive", "m"])
## t = 4.8527, df = 136, p-value = 3.29e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2319669 0.5180282
## sample estimates:
##       cor 
## 0.3841799
cor.test(s2_perf_overall$m, o_min_c_perf_s2)
## 
##  Pearson's product-moment correlation
## 
## data:  s2_perf_overall$m and o_min_c_perf_s2
## t = -0.0051584, df = 136, p-value = 0.9959
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1675348  0.1666748
## sample estimates:
##           cor 
## -0.0004423248
ComparePars(s2_perf_overall$m, o_min_c_perf_s2, use_identity_line = 0)

Invalid relationships to perf diffs

s1_o_inv <- s1_train %>% group_by(ID) %>% summarize(o_inv=mean(overt_invalid)) 
s1_c_inv <- s1_train %>% group_by(ID) %>% summarize(c_inv=mean(cognitive_invalid)) 
assert(all(s1_o_inv$ID==s1_c_inv$ID))
s1_inv_diff <- s1_c_inv$c_inv-s1_o_inv$o_inv
cor.test(o_min_c_perf_s1, s1_inv_diff)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s1 and s1_inv_diff
## t = 0.029401, df = 123, p-value = 0.9766
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1730371  0.1781756
## sample estimates:
##         cor 
## 0.002650977
ComparePars(o_min_c_perf_s1, s1_inv_diff, use_identity_line = 0)

s2_o_inv <- s2_train %>% group_by(ID) %>% summarize(o_inv=mean(prop_invalid_overt)) 
s2_c_inv <- s2_train %>% group_by(ID) %>% summarize(c_inv=mean(prop_invalid_cognitive)) 
assert(all(s2_o_inv$ID==s2_c_inv$ID))
s2_inv_diff <- s2_c_inv$c_inv-s2_o_inv$o_inv
cor.test(o_min_c_perf_s2, s2_inv_diff)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s2 and s2_inv_diff
## t = 2.7091, df = 136, p-value = 0.007615
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.0614962 0.3790481
## sample estimates:
##       cor 
## 0.2262758
ComparePars(s2_inv_diff, o_min_c_perf_s2, use_identity_line = 0, xchar="cog minus overt invalid", ychar="overt minus cog correct")

Percent of pts who actually showing better performance in cognitive

Training and test study 1

# Study 1  
length(which(o_min_c_perf_s1 < 0))/length(o_min_c_perf_s1)
## [1] 0.376
length(which(o_min_c_perf_s1 == 0))/length(o_min_c_perf_s1) # 2.4% show no diff  
## [1] 0.024
length(which(o_min_c_testperf_s1 < 0))/length(o_min_c_testperf_s1)
## [1] 0.328
length(which(o_min_c_testperf_s1 == 0))/length(o_min_c_testperf_s1) 
## [1] 0.104

Training and test study 2

length(which(o_min_c_perf_s2 < 0))/length(o_min_c_perf_s2)
## [1] 0.3913043
length(which(o_min_c_perf_s2 == 0))/length(o_min_c_perf_s2) # < 1% show no diff
## [1] 0.007246377
length(which(o_min_c_testperf_s2 < 0))/length(o_min_c_testperf_s2)
## [1] 0.3550725
length(which(o_min_c_testperf_s2 == 0))/length(o_min_c_testperf_s2) 
## [1] 0.1884058

Modeling Results (Note: model numbers match those from s.R but not those used in paper…)

… because the paper numbering follows the logic of how presented there, and because we culled poor/uninformative/auxiliary models as developing the streamlined final set

Model 1 — Q learner with different decay (phi). Includes eps, CK, block explor parameter

Numbers from before bug fix on the use of prior (3/23/23) dont use: 58986 | 52624 | 57510 | 74680

m1_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior17352.csv")
m1_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior58319.csv")

m1_study1_eb <- rbind(m1_study1_eb_v1, m1_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

m1_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior12123.csv")
m1_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior74336.csv")

m1_study2_eb <- rbind(m1_study2_eb_v1, m1_study2_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

# write.csv(m1_study1_eb,
#           "../model_res/opts_mle_paper_final/best/BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior_merged.csv")
# write.csv(m1_study2_eb,
#           "../model_res/opts_mle_paper_final/best/BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior_merged.csv")
cat("\nBeta Median \n")
## 
## Beta Median
median(m1_study1_eb$beta)
## [1] 2.330957
cat("SD \n")
## SD
sd(m1_study1_eb$beta)
## [1] 5.393476
cat("\nRL LR \n")
## 
## RL LR
median(m1_study1_eb$q_LR)
## [1] 0.5916992
cat("SD \n")
## SD
sd(m1_study1_eb$q_LR)
## [1] 0.2752539
cat("\nphi cog \n")
## 
## phi cog
median(m1_study1_eb$cog_phi)
## [1] 0.1593268
cat("SD \n")
## SD
sd(m1_study1_eb$cog_phi)
## [1] 0.2250586
cat("\nphi overt \n")
## 
## phi overt
median(m1_study1_eb$overt_phi)
## [1] 0.1328623
cat("SD \n")
## SD
sd(m1_study1_eb$overt_phi)
## [1] 0.169769
cat("\nepsilon \n")
## 
## epsilon
median(m1_study1_eb$epsilon)
## [1] 0.005554817
cat("SD \n")
## SD
sd(m1_study1_eb$epsilon)
## [1] 0.02699492
cat("\nChoice LR \n")
## 
## Choice LR
median(m1_study1_eb$choice_LR)
## [1] 0.1198229
cat("SD \n")
## SD
sd(m1_study1_eb$choice_LR)
## [1] 0.1222895
cat("\nES \n")
## 
## ES
median(m1_study1_eb$explor_scalar)
## [1] 0.3783385
cat("SD \n")
## SD
sd(m1_study1_eb$explor_scalar)
## [1] 0.2089571
cat("\nBeta Median \n")
## 
## Beta Median
median(m1_study2_eb$beta)
## [1] 2.692798
cat("SD \n")
## SD
sd(m1_study2_eb$beta)
## [1] 4.910113
cat("\nRL LR \n")
## 
## RL LR
median(m1_study2_eb$q_LR)
## [1] 0.5543408
cat("SD \n")
## SD
sd(m1_study2_eb$q_LR)
## [1] 0.2806
cat("\nphi cog \n")
## 
## phi cog
median(m1_study2_eb$cog_phi)
## [1] 0.2033696
cat("SD \n")
## SD
sd(m1_study2_eb$cog_phi)
## [1] 0.2858507
cat("\nphi overt \n")
## 
## phi overt
median(m1_study2_eb$overt_phi)
## [1] 0.157541
cat("SD \n")
## SD
sd(m1_study2_eb$overt_phi)
## [1] 0.2223507
cat("\nepsilon \n")
## 
## epsilon
median(m1_study2_eb$epsilon)
## [1] 7.844961e-07
cat("SD \n")
## SD
sd(m1_study2_eb$epsilon)
## [1] 0.02668559
cat("\nChoice LR \n")
## 
## Choice LR
median(m1_study2_eb$choice_LR)
## [1] 0.09913484
cat("SD \n")
## SD
sd(m1_study2_eb$choice_LR)
## [1] 0.128331
cat("\nES \n")
## 
## ES
median(m1_study2_eb$explor_scalar)
## [1] 0.2382697
cat("SD \n")
## SD
sd(m1_study2_eb$explor_scalar)
## [1] 0.1951073
ComparePars(m1_study1_eb$cog_phi, m1_study1_eb$overt_phi, 
            "Phi", "Cog", "Overt")

ComparePars(m1_study2_eb$cog_phi, m1_study2_eb$overt_phi, 
            "Phi", "Cog", "Overt")

Relationships at both train and test in studies 1 and 2

assert(s1_perf[s1_perf$condition=="overt", "ID"]==m1_study1_eb$ID)

o_min_c_phi_s1_m1 <- m1_study1_eb$overt_phi - m1_study1_eb$cog_phi

cor.test(o_min_c_perf_s1, o_min_c_phi_s1_m1)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s1 and o_min_c_phi_s1_m1
## t = -6.9629, df = 123, p-value = 1.765e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6469195 -0.3927853
## sample estimates:
##        cor 
## -0.5317171
cor.test(o_min_c_testperf_s1, o_min_c_phi_s1_m1)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s1 and o_min_c_phi_s1_m1
## t = -6.0126, df = 123, p-value = 1.926e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6018399 -0.3284889
## sample estimates:
##        cor 
## -0.4766035
assert(s2_perf[s2_perf$condition=="overt", "ID"]==m1_study2_eb$ID)

o_min_c_phi_s2_m1 <- m1_study2_eb$overt_phi - m1_study2_eb$cog_phi

cor.test(o_min_c_perf_s2, o_min_c_phi_s2_m1)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s2 and o_min_c_phi_s2_m1
## t = -4.7721, df = 136, p-value = 4.642e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5133414 -0.2259169
## sample estimates:
##        cor 
## -0.3787242
cor.test(o_min_c_testperf_s2, o_min_c_phi_s2_m1)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s2 and o_min_c_phi_s2_m1
## t = -7.1401, df = 136, p-value = 5.093e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6339502 -0.3889985
## sample estimates:
##        cor 
## -0.5221611

Descriptives and chi-square tests

Higher median decay in cog in both

# hist(sqrt(m1_study1_eb$cog_phi), breaks=100)
# hist(m1_study1_eb$cog_phi, breaks=100)
# hist(m1_study1_eb$overt_phi, breaks=100)
cat("\n Study 1 cog\n")
## 
##  Study 1 cog
median(m1_study1_eb$cog_phi)
## [1] 0.1593268
cat("\n Study 1 overt\n")
## 
##  Study 1 overt
median(m1_study1_eb$overt_phi)
## [1] 0.1328623
m1_s1_phi <- rbind(
  data.frame("phi"=m1_study1_eb$overt_phi, "cond"="overt"),
  data.frame("phi"=m1_study1_eb$cog_phi, "cond"="cognitive")
)

m1_s2_phi <- rbind(
  data.frame("phi"=m1_study2_eb$overt_phi, "cond"="overt"),
  data.frame("phi"=m1_study2_eb$cog_phi, "cond"="cognitive")
)
m1_s1_phi$Experiment <- "Experiment 1"
m1_s2_phi$Experiment <- "Experiment 2"

m1_phi_ID_comb <- rbind(m1_s1_phi, m1_s2_phi)

phi_medians <- rbind(
  data.frame(m1_s1_phi %>% group_by(cond) %>% summarize(m=median(phi)), "Experiment"="Experiment 1"),
  data.frame(m1_s2_phi %>% group_by(cond) %>% summarize(m=median(phi)), "Experiment"="Experiment 2")
)

Difference higher for means but not as meaningful bc of skew, so plotting median

m1_s1_phi %>% group_by(cond) %>% summarize(m=mean(phi))
## # A tibble: 2 × 2
##   cond          m
##   <chr>     <dbl>
## 1 cognitive 0.235
## 2 overt     0.176
m1_s2_phi %>% group_by(cond) %>% summarize(m=mean(phi))
## # A tibble: 2 × 2
##   cond          m
##   <chr>     <dbl>
## 1 cognitive 0.310
## 2 overt     0.209
phi_plot <-ggplot(phi_medians, aes(x=cond, y=m, color=cond, group=cond)) +
  
  geom_jitter(data = m1_phi_ID_comb, 
              aes(x=cond, y=phi, group=cond), 
              color="black", fill="gray57", size=2, pch=21, width=.15) + 
  geom_bar(stat="identity", fill="white", size=3, alpha=.2) +
  facet_wrap(~ Experiment)  + ga + ap + ft + lp + ylab(TeX('$\\phi') ) + xlab("") + 
  scale_color_manual(values=c("purple", "orange")) + ylim(-.03, 1.1) + tol
phi_plot

#ggsave("../paper/figs/fig_parts/phi_plot.png", phi_plot, height=5, width=10, dpi=700)

Percent higher phi in studies 1 and 2

c_min_o_phi_s1_m1 <- m1_study1_eb$cog_phi - m1_study1_eb$overt_phi
c_min_o_phi_s2_m1 <- m1_study2_eb$cog_phi - m1_study2_eb$overt_phi
table((c_min_o_phi_s1_m1 > 0)*1)[2]/sum(table((c_min_o_phi_s1_m1 > 0)*1))
##     1 
## 0.576
table((c_min_o_phi_s2_m1 > 0)*1)[2]/sum(table((c_min_o_phi_s2_m1 > 0)*1))
##         1 
## 0.6304348
chisq.test(table((c_min_o_phi_s1_m1 > 0)*1))
## 
##  Chi-squared test for given probabilities
## 
## data:  table((c_min_o_phi_s1_m1 > 0) * 1)
## X-squared = 2.888, df = 1, p-value = 0.08924
chisq.test(table((c_min_o_phi_s2_m1 > 0)*1))
## 
##  Chi-squared test for given probabilities
## 
## data:  table((c_min_o_phi_s2_m1 > 0) * 1)
## X-squared = 9.3913, df = 1, p-value = 0.00218
m1_s1_plot_df <- data.frame("train_oc_diff"=o_min_c_perf_s1, "test_oc_diff"=o_min_c_testperf_s1, "phi_co_diff"=c_min_o_phi_s1_m1)
m1_s2_plot_df <- data.frame("train_oc_diff"=o_min_c_perf_s2, "test_oc_diff"=o_min_c_testperf_s2, "phi_co_diff"=c_min_o_phi_s2_m1)

Training plots

a <- ggplot(m1_s1_plot_df, aes(x=phi_co_diff, y=train_oc_diff)) + 
  geom_smooth(method="lm", se=FALSE, size=3, color="black") +
  geom_point(pch=21, fill="gray57", size=6) + 
  stat_cor(method="pearson", size=6, label.y=.5) + 
  ylab("Proportion correct: \nOvert minus cognitive") + 
  xlab(TeX("")) +
  ga + ap + 
  xlim(-.9, 1.05) + ylim(-.7, .65)
a
## `geom_smooth()` using formula = 'y ~ x'

b <- ggplot(m1_s2_plot_df, aes(x=phi_co_diff, y=train_oc_diff)) + 
  geom_smooth(method="lm", se=FALSE, size=3, color="black") +
  geom_point(pch=21, fill="gray57", size=6) + 
  stat_cor(method="pearson", size=6, label.y=.5) + 
  ylab("") + 
  xlab(TeX("")) +
  ga + ap + 
  xlim(-.9, 1.05) + ylim(-.7, .65)
b
## `geom_smooth()` using formula = 'y ~ x'

Test plots

c <- ggplot(m1_s1_plot_df, aes(x=phi_co_diff, y=test_oc_diff)) + 
  geom_smooth(method="lm", se=FALSE, size=3, color="gray57") +
  geom_point(pch=21, fill="gray79", size=6) + 
  stat_cor(method="pearson", size=6, label.y=.5) + 
  ylab("Proportion correct: \nOvert minus cognitive") + 
  xlab(TeX("$\\phi^{Cog}-\\phi^{Overt}")) +
  ga + ap + 
  xlim(-.9, 1.05) + ylim(-.7, .65)

d <- ggplot(m1_s2_plot_df, aes(x=phi_co_diff, y=test_oc_diff)) + 
  geom_smooth(method="lm", se=FALSE, size=3, color="gray57") +
  geom_point(pch=21, fill="gray79", size=6) + 
  stat_cor(method="pearson", size=6, label.y=.5) + 
  ylab("") + 
  xlab(TeX("$\\phi^{Cog}-\\phi^{Overt}")) +
  ga + ap + 
  xlim(-.9, 1.05) + ylim(-.7, .65)
all_vs <- (a + b) / (c + d)
all_vs
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

#ggsave("../paper/figs/fig_parts/phi_vs_propdiff.png", all_vs, height=8, width=14, dpi=700)

Correlation with delay effect

cor.test(m1_study1_eb_v1$cog_phi+m1_study1_eb_v1$overt_phi, ranef(m1_full_delay_cond)$ID$delay_cc1)
## 
##  Pearson's product-moment correlation
## 
## data:  m1_study1_eb_v1$cog_phi + m1_study1_eb_v1$overt_phi and ranef(m1_full_delay_cond)$ID$delay_cc1
## t = -3.7767, df = 123, p-value = 0.0002462
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4712817 -0.1555501
## sample estimates:
##        cor 
## -0.3223523
cor.test(m1_study2_eb_v1$cog_phi+m1_study2_eb_v1$overt_phi, ranef(m2_full_delay_cond)$ID$delay_cc1)
## 
##  Pearson's product-moment correlation
## 
## data:  m1_study2_eb_v1$cog_phi + m1_study2_eb_v1$overt_phi and ranef(m2_full_delay_cond)$ID$delay_cc1
## t = -5.6966, df = 136, p-value = 7.246e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5646058 -0.2933200
## sample estimates:
##        cor 
## -0.4389115
ComparePars((m1_study1_eb_v1$cog_phi+m1_study1_eb_v1$overt_phi), ranef(m1_full_delay_cond)$ID$delay_cc1, use_identity_line = 0)

ComparePars((m1_study2_eb_v1$cog_phi+m1_study2_eb_v1$overt_phi), ranef(m2_full_delay_cond)$ID$delay_cc1, use_identity_line = 0)

For resub, examining whether phi correlates with RT difference

assert(all(s1_tr_diffs_inds$ID==m1_study1_eb$ID))
assert(all(s1_te_diffs_inds$ID==m1_study1_eb$ID))
assert(all(s2_tr_diffs_inds$ID==m1_study2_eb$ID))
assert(all(s2_te_diffs_inds$ID==m1_study2_eb$ID))

ComparePars(m1_study1_eb$cog_phi-m1_study1_eb$overt_phi, s1_tr_diffs_inds$m, use_identity_line = 0)

ComparePars(m1_study1_eb$cog_phi-m1_study1_eb$overt_phi, s1_te_diffs_inds$m, use_identity_line = 0)

ComparePars(m1_study2_eb$cog_phi-m1_study2_eb$overt_phi, s2_tr_diffs_inds$m, use_identity_line = 0)

ComparePars(m1_study2_eb$cog_phi-m1_study2_eb$overt_phi, s2_te_diffs_inds$m, use_identity_line = 0)

cor.test(m1_study1_eb$cog_phi-m1_study1_eb$overt_phi, s1_tr_diffs_inds$m)
## 
##  Pearson's product-moment correlation
## 
## data:  m1_study1_eb$cog_phi - m1_study1_eb$overt_phi and s1_tr_diffs_inds$m
## t = 2.2364, df = 123, p-value = 0.02712
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.02286048 0.36075890
## sample estimates:
##       cor 
## 0.1976744
cor.test(m1_study1_eb$cog_phi-m1_study1_eb$overt_phi, s1_te_diffs_inds$m)
## 
##  Pearson's product-moment correlation
## 
## data:  m1_study1_eb$cog_phi - m1_study1_eb$overt_phi and s1_te_diffs_inds$m
## t = 3.2832, df = 123, p-value = 0.001336
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1139323 0.4376519
## sample estimates:
##       cor 
## 0.2838605
cor.test(m1_study2_eb$cog_phi-m1_study2_eb$overt_phi, s2_tr_diffs_inds$m)
## 
##  Pearson's product-moment correlation
## 
## data:  m1_study2_eb$cog_phi - m1_study2_eb$overt_phi and s2_tr_diffs_inds$m
## t = 4.1166, df = 136, p-value = 6.634e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1755209 0.4736233
## sample estimates:
##       cor 
## 0.3328628
cor.test(m1_study2_eb$cog_phi-m1_study2_eb$overt_phi, s2_te_diffs_inds$m)
## 
##  Pearson's product-moment correlation
## 
## data:  m1_study2_eb$cog_phi - m1_study2_eb$overt_phi and s2_te_diffs_inds$m
## t = 2.6808, df = 136, p-value = 0.008252
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.05914673 0.37702705
## sample estimates:
##       cor 
## 0.2240372

For resub, examining if phi differs based on condition order

unique_IDs <- unique(m1_study1_eb$ID)
for (s in 1:length(unique_IDs)) {
  m1_study1_eb[m1_study1_eb$ID == unique_IDs[s], "first_cond"] <- 
    unique(s1_with_cond_order[s1_with_cond_order$deident_ID == unique_IDs[s], "first_cond"])
}
unique_IDs <- unique(m1_study2_eb$ID)
for (s in 1:length(unique_IDs)) {
  m1_study2_eb[m1_study2_eb$ID == unique_IDs[s], "first_cond"] <- 
    unique(s2_with_cond_order[s2_with_cond_order$de_ID == unique_IDs[s], "first_cond"])
}
m1_study1_eb %>% group_by(first_cond) %>% summarize(om=median(overt_phi), cm=median(cog_phi))
## # A tibble: 2 × 3
##   first_cond       om    cm
##   <chr>         <dbl> <dbl>
## 1 cognTraining  0.118 0.180
## 2 overtTraining 0.148 0.133
m1_study2_eb %>% group_by(first_cond) %>% summarize(om=median(overt_phi), cm=median(cog_phi))
## # A tibble: 2 × 3
##   first_cond    om    cm
##   <chr>      <dbl> <dbl>
## 1 cognitive  0.126 0.354
## 2 overt      0.173 0.141

For resub, show delay distribution

hist(s1_train$delay, breaks=100)

hist(s2_train$delay, breaks=100)

Model validation, main text parts

6/5/23 — reran sims to make sure completely updated/bug fixed version of par results used

s1_train_sim_m1_eb <- 
  rm(sp, 
     "SIM_EMPIRICAL_BAYES_study_1_train_SIT__train_RunMQLearnerDiffDecayToPessPrior57088.csv")

s2_train_sim_m1_eb <- 
  rm(sp, 
    "SIM_EMPIRICAL_BAYES_study_2_train_SIT__train_RunMQLearnerDiffDecayToPessPrior28414.csv")

s1_test_sim_m1_eb <- 
  rm(sp, 
     "SIM_EMPIRICAL_BAYES_study_1_train_SIT__sit_RunMQLearnerDiffDecayToPessPrior57088.csv")

s2_test_sim_m1_eb <- 
  rm(sp, 
    "SIM_EMPIRICAL_BAYES_study_2_train_SIT__sit_RunMQLearnerDiffDecayToPessPrior28414.csv")
a <- AltOneSimEmpPlotRep(AltPlotSimTrainingCurves(emp_df=s1_train, s1_train_sim_m1_eb))
b <- AltOneSimEmpPlotRep(AltPlotSimTrainingCurves(emp_df=s2_train, s2_train_sim_m1_eb))
c <- PlotSimEmpTest(emp_df=s1_sit, sim_df=s1_test_sim_m1_eb)
d <- PlotSimEmpTest(emp_df=s2_sit, sim_df=s2_test_sim_m1_eb)
a

b

sim_test_comb <- c + d 
#(temporarily flipped tol -> lp to get the legend for paper)  
sim_test_comb

#ggsave("../paper/figs/fig_parts/train_sim_experiment1_plot.png", a, height=7, width=12, dpi=700)
ggsave("../paper/figs/fig_parts/test_sim_experiment1_plot.png", c, height=5, width=11, dpi=700)
ggsave("../paper/figs/fig_parts/test_sim_experiment2_plot.png", d, height=5, width=11, dpi=700)

For supplemental because of space constraints

#ggsave("../paper/figs/fig_parts/train_sim_experiment2_plot.png", b, height=7, width=12, dpi=700)

Model 3 — Same as m1 but not letting decay vary

Prior to vactoin — bug fixed files swapped in (see to_do notes)

m3_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior35433.csv")
m3_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior22212.csv")

m3_study1_eb <- rbind(m3_study1_eb_v1, m3_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

m3_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior29818.csv")
m3_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior56742.csv")

m3_study2_eb <- rbind(m3_study2_eb_v1,m3_study2_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

# write.csv(m3_study1_eb,
#           "../model_res/opts_mle_paper_final/best/BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior_merged.csv")
# write.csv(m3_study2_eb,
#           "../model_res/opts_mle_paper_final/best/BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior_merged.csv")

Worse in AIC

ComparePars(m3_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m3", "m1")

ComparePars(m3_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m3", "m1")

sum(m3_study1_eb$AIC-m1_study1_eb$AIC)
## [1] 134.026
sum(m3_study2_eb$AIC-m1_study2_eb$AIC)
## [1] 147.2234

Sanity check that the amount of model improvement when letting decay varies correlates with the diff in phi between conditions

assert(m1_study1_eb$ID==m3_study1_eb$ID)
assert(m1_study2_eb$ID==m3_study2_eb$ID)
cor.test(m3_study1_eb$AIC-m1_study1_eb$AIC, abs(m1_study1_eb$overt_phi - m1_study1_eb$cog_phi), method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  m3_study1_eb$AIC - m1_study1_eb$AIC and abs(m1_study1_eb$overt_phi - m1_study1_eb$cog_phi)
## S = 107038, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.6711582
cor.test(m3_study2_eb$AIC-m1_study2_eb$AIC, abs(m1_study2_eb$overt_phi - m1_study2_eb$cog_phi), method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  m3_study2_eb$AIC - m1_study2_eb$AIC and abs(m1_study2_eb$overt_phi - m1_study2_eb$cog_phi)
## S = 163526, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.6266436

(Function prints Pearson’s r but more appropriate Spearman’s \(\rho\) reported one cell above)

ComparePars(m1_study1_eb$AIC-m3_study1_eb$AIC, abs(m1_study1_eb$overt_phi - m1_study1_eb$cog_phi), "", "AIC change", "Overt min Cog phi", use_identity_line = 0)

ComparePars(m1_study2_eb$AIC-m3_study2_eb$AIC, abs(m1_study2_eb$overt_phi - m1_study2_eb$cog_phi), "", "AIC change", "Overt min Cog phi", use_identity_line = 0)

Model 4 — Same as m3 but to 0 inits

m4_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits47402.csv")
m4_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits84391.csv")

m4_study1_eb <- rbind(m4_study1_eb_v1, m4_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

m4_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits41957.csv")
m4_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits37660.csv")

m4_study2_eb <- rbind(m4_study2_eb_v1, 
                      m4_study2_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

# write.csv(m4_study1_eb,
#           "../model_res/opts_mle_paper_final/best/BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits_merged.csv")
# write.csv(m4_study2_eb,
#           "../model_res/opts_mle_paper_final/best/BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits_merged.csv")

Worse in AIC than pessimistic priors

ComparePars(m4_study1_eb$AIC, m3_study1_eb$AIC, "AIC", "m3", "m1")

ComparePars(m4_study2_eb$AIC, m3_study2_eb$AIC, "AIC", "m3", "m1")

sum(m4_study1_eb$AIC-m3_study1_eb$AIC)
## [1] 120.8681
sum(m4_study2_eb$AIC-m3_study2_eb$AIC)
## [1] 446.3359

Model 5 — Same as m1 but lets beta vary instead of decay

4/6 — bug fixed

m5_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffBeta67330.csv")
m5_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffBeta37484.csv")

m5_study1_eb <- rbind(m5_study1_eb_v1, m5_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

m5_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffBeta76602.csv")
m5_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffBeta11129.csv")

m5_study2_eb <- rbind(m5_study2_eb_v1, m5_study2_eb_v2) %>% group_by(ID) %>% slice(which.min(nll))

Worse in AIC than diff decay

ComparePars(m5_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m5", "m1")

ComparePars(m5_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m5", "m1")

sum(m5_study1_eb$AIC-m1_study1_eb$AIC)
## [1] 672.2797
sum(m5_study2_eb$AIC-m1_study2_eb$AIC)
## [1] 857.6129

Model 6 — Same as m1 but lets learning rate vary instead of decay

4/6 — bug fix files swapped in

m6_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffLR12095.csv")
m6_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffLR13604.csv")

m6_study1_eb <- rbind(m6_study1_eb_v1, m6_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

m6_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffLR69123.csv")
m6_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffLR43205.csv")

m6_study2_eb <- rbind(m6_study2_eb_v1, m6_study2_eb_v2) %>% group_by(ID) %>% slice(which.min(nll))

Worse in AIC than decay allowed to vary by condition

ComparePars(m6_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m6", "m1")

ComparePars(m6_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m6", "m1")

sum(m6_study1_eb$AIC-m1_study1_eb$AIC)
## [1] 131.9173
sum(m6_study2_eb$AIC-m1_study2_eb$AIC)
## [1] 173.9311

Higher overt learning rate in both

median(m6_study1_eb$cog_q_LR)
## [1] 0.5048065
median(m6_study1_eb$overt_q_LR)
## [1] 0.556566
median(m6_study2_eb$cog_q_LR)
## [1] 0.4577678
median(m6_study1_eb$overt_q_LR)
## [1] 0.556566

Learning rate correlated across conditions

ComparePars(m6_study1_eb$cog_q_LR, m6_study1_eb$overt_q_LR, 
            "q_LR", "Cog", "Overt")

ComparePars(m6_study2_eb$cog_q_LR, m6_study2_eb$overt_q_LR, 
            "q_LR", "Cog", "Overt")

Correlated with perf diffs but not as strongly as phi

assert(s1_perf[s1_perf$condition=="overt", "ID"]==m6_study1_eb$ID)

o_min_c_q_LR_s1_m6 <- m6_study1_eb$overt_q_LR - m6_study1_eb$cog_q_LR

cor.test(o_min_c_perf_s1, o_min_c_q_LR_s1_m6)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s1 and o_min_c_q_LR_s1_m6
## t = 4.361, df = 123, p-value = 2.708e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2034049 0.5088491
## sample estimates:
##       cor 
## 0.3659413
cor.test(o_min_c_testperf_s1, o_min_c_q_LR_s1_m6)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s1 and o_min_c_q_LR_s1_m6
## t = 2.6939, df = 123, p-value = 0.008049
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.06303996 0.39525878
## sample estimates:
##       cor 
## 0.2360345
assert(s2_perf[s2_perf$condition=="overt", "ID"]==m6_study2_eb$ID)

o_min_c_q_LR_s2_m6 <- m6_study2_eb$overt_q_LR - m6_study2_eb$cog_q_LR

cor.test(o_min_c_perf_s2, o_min_c_q_LR_s2_m6)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s2 and o_min_c_q_LR_s2_m6
## t = 3.83, df = 136, p-value = 0.0001949
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1528938 0.4553870
## sample estimates:
##       cor 
## 0.3120266
cor.test(o_min_c_testperf_s2, o_min_c_q_LR_s2_m6)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s2 and o_min_c_q_LR_s2_m6
## t = 5.4029, df = 136, p-value = 2.849e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2724024 0.5489175
## sample estimates:
##      cor 
## 0.420372
chisq.test(table((o_min_c_q_LR_s1_m6 > 0)*1))
## 
##  Chi-squared test for given probabilities
## 
## data:  table((o_min_c_q_LR_s1_m6 > 0) * 1)
## X-squared = 6.728, df = 1, p-value = 0.009491
chisq.test(table((o_min_c_q_LR_s2_m6 > 0)*1))
## 
##  Chi-squared test for given probabilities
## 
## data:  table((o_min_c_q_LR_s2_m6 > 0) * 1)
## X-squared = 1.8551, df = 1, p-value = 0.1732

Model 11 — Simple q-learner model

m11_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearner24681.csv")
m11_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearner77403.csv")

m11_study1_eb <- rbind(m11_study1_eb_v1,m11_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))


m11_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearner18790.csv")
m11_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearner26119.csv")

m11_study2_eb <- rbind(m11_study2_eb_v1, m11_study2_eb_v2) %>% group_by(ID) %>% slice(which.min(nll))

Much worse in AIC to not include decay

ComparePars(m11_study1_eb$AIC, m3_study1_eb$AIC, "AIC", "m11", "m3")

ComparePars(m11_study2_eb$AIC, m3_study2_eb$AIC, "AIC", "m11", "m3")

sum(m11_study1_eb$AIC-m3_study1_eb$AIC)
## [1] 918.2861
sum(m11_study2_eb$AIC-m3_study2_eb$AIC)
## [1] 1225.945

Model 12 — Simple Bayes model

m12_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMBayesLearner13415.csv")
m12_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMBayesLearner42048.csv")

m12_study1_eb <- rbind(m12_study1_eb_v1,m12_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

m12_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMBayesLearner22977.csv")
m12_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMBayesLearner23536.csv")

m12_study2_eb <- rbind(m12_study2_eb_v1, m12_study2_eb_v2) %>% group_by(ID) %>% slice(which.min(nll))

In turn much worse than the simple Q-learning model

ComparePars(m12_study1_eb$AIC, m11_study1_eb$AIC, "AIC", "m12", "m11")

ComparePars(m12_study2_eb$AIC, m11_study2_eb$AIC, "AIC", "m12", "m11")

sum(m12_study1_eb$AIC-m11_study1_eb$AIC)
## [1] 622.5435
sum(m12_study2_eb$AIC-m11_study2_eb$AIC)
## [1] 1012.796

Model 13 — M12 but allowing beta to vary

m13_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMBayesLearnerDiffBeta60557.csv")
m13_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMBayesLearnerDiffBeta85855.csv")

m13_study1_eb <- rbind(m13_study1_eb_v1, m13_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

m13_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMBayesLearnerDiffBeta29859.csv")
m13_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMBayesLearnerDiffBeta29859.csv")

m13_study2_eb <- rbind(m13_study2_eb_v1, m13_study2_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

Way worse than m1

ComparePars(m13_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m13", "m1")

ComparePars(m13_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m13", "m1")

sum(m13_study1_eb$AIC-m1_study1_eb$AIC)
## [1] 2081.633
sum(m13_study2_eb$AIC-m1_study2_eb$AIC)
## [1] 2997.104

Higher median beta in overt

median(m13_study1_eb$cog_beta)
## [1] 1.790819
median(m13_study1_eb$overt_beta)
## [1] 1.876175
median(m13_study2_eb$cog_beta)
## [1] 1.893116
median(m13_study2_eb$overt_beta)
## [1] 2.099269

Modest correlations between beta and perf

assert(s1_perf[s1_perf$condition=="overt", "ID"]==m13_study1_eb$ID)

o_min_c_beta_s1_m13 <- m13_study1_eb$overt_beta - m13_study1_eb$cog_beta

cor.test(o_min_c_perf_s1, o_min_c_beta_s1_m13)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s1 and o_min_c_beta_s1_m13
## t = 3.5966, df = 123, p-value = 0.0004653
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1404791 0.4592086
## sample estimates:
##       cor 
## 0.3084768
cor.test(o_min_c_testperf_s1, o_min_c_beta_s1_m13)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s1 and o_min_c_beta_s1_m13
## t = 2.7331, df = 123, p-value = 0.007199
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.06646281 0.39815509
## sample estimates:
##       cor 
## 0.2392776
assert(s2_perf[s2_perf$condition=="overt", "ID"]==m13_study2_eb$ID)

o_min_c_beta_s2_m13 <- m13_study2_eb$overt_beta - m13_study2_eb$cog_beta

cor.test(o_min_c_perf_s2, o_min_c_beta_s2_m13)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s2 and o_min_c_beta_s2_m13
## t = 2.8955, df = 136, p-value = 0.004412
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.07696563 0.39227998
## sample estimates:
##       cor 
## 0.2409713
cor.test(o_min_c_testperf_s2, o_min_c_beta_s2_m13)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s2 and o_min_c_beta_s2_m13
## t = 1.9265, df = 136, p-value = 0.05613
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.004233115  0.321339865
## sample estimates:
##      cor 
## 0.162987

Model 14 — Hybrid model — integrative plus bayes — bayes contribution varies by cond

Alt refers to this being the simpler version with uncertainty factored in at the action (arm) level rather than a weighted uncertainty computed — this and the more complex version were v similar in terms of fit

m14_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayesAlt39650.csv")
m14_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayesAlt71571.csv")

m14_study1_eb <- rbind(m14_study1_eb_v1, m14_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

m14_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayesAlt38226.csv")
m14_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayesAlt28286.csv")

m14_study2_eb <- rbind(m14_study2_eb_v1, m14_study2_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

Original version with the more complicated weighting scheme

m14_study1_eb_v1_old <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayes25885.csv")
m14_study1_eb_v2_old <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayes65430.csv")

m14_study1_eb_old <- rbind(m14_study1_eb_v1_old, m14_study1_eb_v2_old) %>% 
  group_by(ID) %>% slice(which.min(nll))

m14_study2_eb_v1_old <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayes30853.csv")
m14_study2_eb_v2_old <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayes58630.csv")

m14_study2_eb_old <- rbind(m14_study2_eb_v1_old, m14_study2_eb_v2_old) %>% 
  group_by(ID) %>% slice(which.min(nll))

As expected, these are very similar, with the new version a tad better

ComparePars(m14_study1_eb$AIC, m14_study1_eb_old$AIC, "AIC", "m14", "m14 old")

sum(m14_study1_eb$AIC-m14_study1_eb_old$AIC)
## [1] -18.67696
ComparePars(m14_study2_eb$AIC, m14_study2_eb_old$AIC, "AIC", "m14", "m14 old")

sum(m14_study2_eb$AIC-m14_study2_eb_old$AIC)
## [1] -43.19563

Worse AIC (and nll) compared to m1 (not a nested model bc doesn’t have diff decay)

ComparePars(m14_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m14", "m1")

ComparePars(m14_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m14", "m1")

sum(m14_study1_eb$AIC-m1_study1_eb$AIC)
## [1] 448.2412
sum(m14_study2_eb$AIC-m1_study2_eb$AIC)
## [1] 486.6508
sum(m14_study1_eb$nll-m1_study1_eb$nll)
## [1] 99.12058
sum(m14_study2_eb$nll-m1_study2_eb$nll)
## [1] 105.3254

Higher median rho in overt

median(m14_study1_eb$rho_cog)
## [1] 0.1651074
median(m14_study1_eb$rho_overt)
## [1] 0.2611568
median(m14_study2_eb$rho_cog)
## [1] 0.1066345
median(m14_study2_eb$rho_overt)
## [1] 0.2367347

Correlations w perf but not as strongly as phi

assert(s1_perf[s1_perf$condition=="overt", "ID"]==m14_study1_eb$ID)

o_min_c_beta_s1_m14 <- m14_study1_eb$rho_overt - m14_study1_eb$rho_cog

cor.test(o_min_c_perf_s1, o_min_c_beta_s1_m14)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s1 and o_min_c_beta_s1_m14
## t = 4.2195, df = 123, p-value = 4.71e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1919767 0.4999821
## sample estimates:
##       cor 
## 0.3555962
cor.test(o_min_c_testperf_s1, o_min_c_beta_s1_m14)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s1 and o_min_c_beta_s1_m14
## t = 3.3007, df = 123, p-value = 0.001262
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1154251 0.4388741
## sample estimates:
##       cor 
## 0.2852507
assert(s2_perf[s2_perf$condition=="overt", "ID"]==m14_study2_eb$ID)

o_min_c_beta_s2_m14 <- m14_study2_eb$rho_overt - m14_study2_eb$rho_cog

cor.test(o_min_c_perf_s2, o_min_c_beta_s2_m14)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s2 and o_min_c_beta_s2_m14
## t = 5.1757, df = 136, p-value = 7.977e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2558970 0.5363993
## sample estimates:
##       cor 
## 0.4056554
cor.test(o_min_c_testperf_s2, o_min_c_beta_s2_m14)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s2 and o_min_c_beta_s2_m14
## t = 3.4613, df = 136, p-value = 0.0007187
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1232954 0.4311429
## sample estimates:
##       cor 
## 0.2845379

Model 15 — Hybrid model — integrative plus bayes — decay contribution varies by cond

Alternate version with the simpler action-wise weighting (see m14 description)

m15_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayesAlt33317.csv")
m15_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayesAlt35801.csv")

m15_study1_eb <- rbind(m15_study1_eb_v1, m15_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))


m15_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayesAlt64665.csv")
m15_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayesAlt54458.csv")

m15_study2_eb <- rbind(m15_study2_eb_v1, m15_study2_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

Original version with the more complicated weighting scheme

bug numbers: 75340 81994 63931 40889

m15_study1_eb_v1_old <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayes50404.csv")
m15_study1_eb_v2_old <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayes56900.csv")

m15_study1_eb_old <- rbind(m15_study1_eb_v1_old, m15_study1_eb_v2_old) %>% 
  group_by(ID) %>% slice(which.min(nll))


m15_study2_eb_v1_old <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayes68726.csv")
m15_study2_eb_v2_old <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayes30103.csv")

m15_study2_eb_old <- rbind(m15_study2_eb_v1_old, m15_study2_eb_v2_old) %>% 
  group_by(ID) %>% slice(which.min(nll))

As expected, again very similar — so using new/simpler one going forward

ComparePars(m15_study1_eb$AIC, m15_study1_eb_old$AIC, "AIC", "m15", "m15 old")

sum(m15_study1_eb$AIC-m15_study1_eb_old$AIC)
## [1] 11.82675
ComparePars(m15_study2_eb$AIC, m15_study2_eb_old$AIC, "AIC", "m15", "m15 old")

sum(m15_study2_eb$AIC-m15_study2_eb_old$AIC)
## [1] -53.60774

Worse AIC compared to m1

ComparePars(m15_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m15", "m1")

ComparePars(m15_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m15", "m1")

sum(m15_study1_eb$AIC-m1_study1_eb$AIC)
## [1] 169.3991
sum(m15_study2_eb$AIC-m1_study2_eb$AIC)
## [1] 127.4026
# sum(m15_study1_eb$nll-m1_study1_eb$nll)
# sum(m15_study2_eb$nll-m1_study2_eb$nll)

Higher median decay in cog in both

median(m15_study1_eb$cog_phi)
## [1] 0.2204569
median(m15_study1_eb$overt_phi)
## [1] 0.1742236
median(m15_study2_eb$cog_phi)
## [1] 0.3308409
median(m15_study1_eb$overt_phi)
## [1] 0.1742236
ComparePars(m15_study1_eb$cog_phi, m15_study1_eb$overt_phi, 
            "Phi", "Cog", "Overt")

ComparePars(m15_study2_eb$cog_phi, m15_study2_eb$overt_phi, 
            "Phi", "Cog", "Overt")

Similar perf correlations as m1

assert(s1_perf[s1_perf$condition=="overt", "ID"]==m15_study1_eb$ID)

o_min_c_phi_s1_m15 <- m15_study1_eb$overt_phi - m15_study1_eb$cog_phi

cor.test(o_min_c_perf_s1, o_min_c_phi_s1_m15)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s1 and o_min_c_phi_s1_m15
## t = -7.1705, df = 123, p-value = 6.097e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6560039 -0.4060503
## sample estimates:
##       cor 
## -0.542943
cor.test(o_min_c_testperf_s1, o_min_c_phi_s1_m15)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s1 and o_min_c_phi_s1_m15
## t = -6.1734, df = 123, p-value = 8.902e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6098839 -0.3397791
## sample estimates:
##        cor 
## -0.4863663
assert(s2_perf[s2_perf$condition=="overt", "ID"]==m15_study2_eb$ID)

o_min_c_phi_s2_m15 <- m15_study2_eb$overt_phi - m15_study2_eb$cog_phi

cor.test(o_min_c_perf_s2, o_min_c_phi_s2_m15)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s2 and o_min_c_phi_s2_m15
## t = -3.8478, df = 136, p-value = 0.0001825
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4565342 -0.1543079
## sample estimates:
##        cor 
## -0.3133333
cor.test(o_min_c_testperf_s2, o_min_c_phi_s2_m15)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s2 and o_min_c_phi_s2_m15
## t = -6.5236, df = 136, p-value = 1.247e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6058791 -0.3496199
## sample estimates:
##        cor 
## -0.4882024

Model 27 — Fixing both ES and epsilon at the group level

m27_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed46645.csv")
m27_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed33416.csv")

m27_study1_eb <- rbind(m27_study1_eb_v1, m27_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

m27_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed20328.csv")
m27_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed44601.csv")

m27_study2_eb <- rbind(m27_study2_eb_v1, m27_study2_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

# write.csv(m27_study1_eb,
#           "../model_res/opts_mle_paper_final/best/BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed_merged.csv")
# write.csv(m27_study2_eb,
#           "../model_res/opts_mle_paper_final/best/BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed_merged.csv")

Does worsen NLL substantially

ComparePars(m27_study1_eb$nll, m1_study1_eb$nll, "nll", "m27", "m1")

ComparePars(m27_study2_eb$nll, m1_study2_eb$nll, "nll", "m27", "m1")

sum(m27_study1_eb$nll-m1_study1_eb$nll)
## [1] 259.0197
sum(m27_study2_eb$nll-m1_study2_eb$nll)
## [1] 367.588

Still higher median decay in cog in both

median(m27_study1_eb$cog_phi)
## [1] 0.1637271
median(m27_study1_eb$overt_phi)
## [1] 0.1337206
median(m27_study2_eb$cog_phi)
## [1] 0.2304866
median(m27_study1_eb$overt_phi)
## [1] 0.1337206

Strong relationships at both train and test in studies 1 and 2

assert(s1_perf[s1_perf$condition=="overt", "ID"]==m27_study1_eb$ID)

o_min_c_phi_s1_m27 <- m27_study1_eb$overt_phi - m27_study1_eb$cog_phi

cor.test(o_min_c_perf_s1, o_min_c_phi_s1_m27)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s1 and o_min_c_phi_s1_m27
## t = -6.9726, df = 123, p-value = 1.68e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6473518 -0.3934142
## sample estimates:
##        cor 
## -0.5322504
cor.test(o_min_c_testperf_s1, o_min_c_phi_s1_m27)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s1 and o_min_c_phi_s1_m27
## t = -6.4826, df = 123, p-value = 1.966e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6248649 -0.3610156
## sample estimates:
##       cor 
## -0.504631
assert(s2_perf[s2_perf$condition=="overt", "ID"]==m27_study2_eb$ID)

o_min_c_phi_s2_m27 <- m27_study2_eb$overt_phi - m27_study2_eb$cog_phi

cor.test(o_min_c_perf_s2, o_min_c_phi_s2_m27)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_perf_s2 and o_min_c_phi_s2_m27
## t = -4.8922, df = 136, p-value = 2.775e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5203133 -0.2349247
## sample estimates:
##        cor 
## -0.3868433
cor.test(o_min_c_testperf_s2, o_min_c_phi_s2_m27)
## 
##  Pearson's product-moment correlation
## 
## data:  o_min_c_testperf_s2 and o_min_c_phi_s2_m27
## t = -7.114, df = 136, p-value = 5.847e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6328069 -0.3873771
## sample estimates:
##       cor 
## -0.520771

59 and 62% of pts show a higher phi in overt

table((o_min_c_phi_s1_m27 > 0)*1)
## 
##  0  1 
## 73 52
table((o_min_c_phi_s2_m27 > 0)*1)
## 
##  0  1 
## 86 52
table((o_min_c_phi_s1_m27 > 0))[1]/sum(table((o_min_c_phi_s1_m27 > 0)))
## FALSE 
## 0.584
table((o_min_c_phi_s2_m27 > 0))[1]/sum(table((o_min_c_phi_s2_m27 > 0)))
##     FALSE 
## 0.6231884
chisq.test(table(if_else((o_min_c_phi_s1_m27 > 0) == TRUE, 1, 0)))
## 
##  Chi-squared test for given probabilities
## 
## data:  table(if_else((o_min_c_phi_s1_m27 > 0) == TRUE, 1, 0))
## X-squared = 3.528, df = 1, p-value = 0.06034
chisq.test(table(if_else((o_min_c_phi_s2_m27 > 0) == TRUE, 1, 0)))
## 
##  Chi-squared test for given probabilities
## 
## data:  table(if_else((o_min_c_phi_s2_m27 > 0) == TRUE, 1, 0))
## X-squared = 8.3768, df = 1, p-value = 0.0038

Model 29 — Model 1 but with no choice kernel

m29_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK21030.csv")
m29_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK89190.csv")
# Replace when back  
#m29_study1_eb <- m29_study1_eb_v1
m29_study1_eb <- rbind(m29_study1_eb_v1, m29_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

m29_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK64103.csv")
m29_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK74167.csv")

# Delete when final one is back  
m29_study2_eb <- m29_study2_eb_v1
  
m29_study2_eb <- rbind(m29_study2_eb_v1, m29_study2_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

# write.csv(m29_study1_eb,
#           "../model_res/opts_mle_paper_final/best/BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK_merged.csv")
# write.csv(m29_study2_eb,
#           "../model_res/opts_mle_paper_final/best/BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK_merged.csv")

Dramatically worse fit

ComparePars(m29_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m29", "m1")

ComparePars(m29_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m29", "m1")

sum(m29_study1_eb$AIC-m1_study1_eb$AIC)
## [1] 1815.325
sum(m29_study2_eb$AIC-m1_study2_eb$AIC)
## [1] 2777.038

Model Comparison

AIC_diff_df <-
 rbind(
    data.frame("diff"=m12_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Bayes"),  
    data.frame("diff"=m12_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Bayes"),
    
    data.frame("diff"=m4_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decayto0"),  
    data.frame("diff"=m4_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decayto0"),  
  
    data.frame("diff"=m3_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decay"),  
    data.frame("diff"=m3_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decay"), 
    
    data.frame("diff"=m5_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decay_diff_beta"),  
    data.frame("diff"=m5_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decay_diff_beta"),
    
    data.frame("diff"=m13_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Bayes_diff_beta"),  
    data.frame("diff"=m13_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Bayes_diff_beta"),
    
    data.frame("diff"=m6_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decay_diff_lr"),  
    data.frame("diff"=m6_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decay_diff_lr"),
    
    data.frame("diff"=m14_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decay_diff_bayes"),  
    data.frame("diff"=m14_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decay_diff_bayes"),
    
    data.frame("diff"=m15_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_diff_decay_bayes"),  
    data.frame("diff"=m15_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_diff_decay_bayes"),

    data.frame("diff"=m1_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decay_diff_phi"),  
    data.frame("diff"=m1_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decay_diff_phi")
 )

AIC_diff_df$study <- as.factor(AIC_diff_df$study)
AIC_diff_summ <- AIC_diff_df %>% group_by(model, study) %>% summarize(m=mean(diff))
## `summarise()` has grouped output by 'model'. You can override using the
## `.groups` argument.
sum(m1_study1_eb$AIC-m3_study1_eb$AIC)
## [1] -134.026
sum(m1_study2_eb$AIC-m3_study2_eb$AIC)
## [1] -147.2234

Confirm min in both studies is diff phi

AIC_diff_summ[AIC_diff_summ$study == 1, ] %>% arrange(m)
## # A tibble: 9 × 3
## # Groups:   model [9]
##   model              study     m
##   <chr>              <fct> <dbl>
## 1 Q_decay_diff_phi   1     -8.42
## 2 Q_decay_diff_lr    1     -7.36
## 3 Q_decay            1     -7.35
## 4 Q_diff_decay_bayes 1     -7.06
## 5 Q_decayto0         1     -6.38
## 6 Q_decay_diff_bayes 1     -4.83
## 7 Q_decay_diff_beta  1     -3.04
## 8 Bayes              1      4.98
## 9 Bayes_diff_beta    1      8.23
AIC_diff_summ[AIC_diff_summ$study == 2, ] %>% arrange(m)
## # A tibble: 9 × 3
## # Groups:   model [9]
##   model              study     m
##   <chr>              <fct> <dbl>
## 1 Q_decay_diff_phi   2     -9.95
## 2 Q_diff_decay_bayes 2     -9.03
## 3 Q_decay            2     -8.88
## 4 Q_decay_diff_lr    2     -8.69
## 5 Q_decay_diff_bayes 2     -6.42
## 6 Q_decayto0         2     -5.65
## 7 Q_decay_diff_beta  2     -3.74
## 8 Bayes              2      7.34
## 9 Bayes_diff_beta    2     11.8
levels <-
  c(
    "Bayes_diff_beta",
    "Bayes", 
    "Q_decayto0", 
    "Q_decay", 
    "Q_decay_diff_beta",  
    "Q_decay_diff_lr", 
    "Q_decay_diff_bayes", 
    "Q_diff_decay_bayes", 
    "Q_decay_diff_phi"
  )

assert(length(levels)==9)

AIC_diff_summ$model <- factor(AIC_diff_summ$model, levels = levels)
AIC_diff_df$model <- factor(AIC_diff_df$model, levels = levels)
mp1 <- ggplot(AIC_diff_summ %>% filter(study == 1), aes(x=model, y=m, group=model, fill=model)) + 
  geom_hline(yintercept = 0, size=4) +
  geom_hline(yintercept = seq(-50, 60, 10), color="gray57", alpha=.8, size=.65) +
  geom_jitter(data=AIC_diff_df %>% filter(study == 1), aes(x=model, y=diff, group=model), 
                fill="gray57", pch=21, width=.3, size=4, alpha=.08) + 
  geom_bar(stat="identity", color="black", alpha=.8, size=2) + 
  ylab(TeX('$\\Delta\ AIC$ from simple Q-learner')) +
  tol + ga + ap + tp +
  ggtitle("Study 1") +
  #facet_wrap(~ study) +
  xlab("") +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())  +
    theme(axis.title.y=element_text(size=22)) + ylim(-55, 65) + coord_flip()
mp2 <- ggplot(AIC_diff_summ %>% filter(study == 2), aes(x=model, y=m, group=model, fill=model)) + 
  geom_hline(yintercept = 0, size=4) +
  geom_hline(yintercept = seq(-50, 60, 10), color="gray57", alpha=.8, size=.65) +
  geom_jitter(data=AIC_diff_df %>% filter(study == 2), aes(x=model, y=diff, group=model), 
                fill="gray57", pch=21, width=.3, size=4, alpha=.08) + 
  geom_bar(stat="identity", color="black", alpha=.8, size=2) + 
  ylab(TeX('$\\Delta\ AIC$ from simple Q-learner')) +
  tol + ga + ap + tp +
  ggtitle("Study 2") +
  xlab("") +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())  +
  theme(axis.title.y=element_text(size=22)) + ylim(-55, 65) + coord_flip()
combined_aic <- mp1 + mp2
combined_aic

#ggsave("../paper/figs/fig_parts/aic_plot.png", combined_aic, height=5, width=10, dpi=700)
AIC_diff_df %>% group_by(model, study) %>% summarize(m=mean(diff)) %>% arrange(m)
## `summarise()` has grouped output by 'model'. You can override using the
## `.groups` argument.
## # A tibble: 18 × 3
## # Groups:   model [9]
##    model              study     m
##    <fct>              <fct> <dbl>
##  1 Q_decay_diff_phi   2     -9.95
##  2 Q_diff_decay_bayes 2     -9.03
##  3 Q_decay            2     -8.88
##  4 Q_decay_diff_lr    2     -8.69
##  5 Q_decay_diff_phi   1     -8.42
##  6 Q_decay_diff_lr    1     -7.36
##  7 Q_decay            1     -7.35
##  8 Q_diff_decay_bayes 1     -7.06
##  9 Q_decay_diff_bayes 2     -6.42
## 10 Q_decayto0         1     -6.38
## 11 Q_decayto0         2     -5.65
## 12 Q_decay_diff_bayes 1     -4.83
## 13 Q_decay_diff_beta  2     -3.74
## 14 Q_decay_diff_beta  1     -3.04
## 15 Bayes              1      4.98
## 16 Bayes              2      7.34
## 17 Bayes_diff_beta    1      8.23
## 18 Bayes_diff_beta    2     11.8

Get AICs for supplemental table

cat(round(sum(m1_study1_eb$AIC), 0), "\n")
## 42125
cat(round(sum(m1_study2_eb$AIC), 0))
## 42798
cat(round(sum(m3_study1_eb$AIC), 0), "\n")
## 42259
cat(round(sum(m3_study2_eb$AIC), 0))
## 42945
cat(round(sum(m4_study1_eb$AIC), 0), "\n")
## 42380
cat(round(sum(m4_study2_eb$AIC), 0))
## 43391
cat(round(sum(m5_study1_eb$AIC), 0), "\n")
## 42797
cat(round(sum(m5_study2_eb$AIC), 0))
## 43655
cat(round(sum(m6_study1_eb$AIC), 0), "\n")
## 42257
cat(round(sum(m6_study2_eb$AIC), 0))
## 42972
cat(round(sum(m13_study1_eb$AIC), 0), "\n")
## 44206
cat(round(sum(m13_study2_eb$AIC), 0))
## 45795
cat(round(sum(m14_study1_eb$AIC), 0), "\n")
## 42573
cat(round(sum(m14_study2_eb$AIC), 0))
## 43284
cat(round(sum(m15_study1_eb$AIC), 0), "\n")
## 42294
cat(round(sum(m15_study2_eb$AIC), 0))
## 42925
cat(round(sum(m11_study1_eb$AIC), 0), "\n")
## 43177
cat(round(sum(m11_study2_eb$AIC), 0))
## 44171
cat(round(sum(m12_study1_eb$AIC), 0), "\n")
## 43800
cat(round(sum(m12_study2_eb$AIC), 0))
## 45184